1 /* IEEE754 floating point arithmetic 2 * single precision 3 */ 4 /* 5 * MIPS floating point support 6 * Copyright (C) 1994-2000 Algorithmics Ltd. 7 * http://www.algor.co.uk 8 * 9 * ######################################################################## 10 * 11 * This program is free software; you can distribute it and/or modify it 12 * under the terms of the GNU General Public License (Version 2) as 13 * published by the Free Software Foundation. 14 * 15 * This program is distributed in the hope it will be useful, but WITHOUT 16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 18 * for more details. 19 * 20 * You should have received a copy of the GNU General Public License along 21 * with this program; if not, write to the Free Software Foundation, Inc., 22 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. 23 * 24 * ######################################################################## 25 */ 26 27 28 #include "ieee754sp.h" 29 30 int ieee754sp_class(ieee754sp x) 31 { 32 COMPXSP; 33 EXPLODEXSP; 34 return xc; 35 } 36 37 int ieee754sp_isnan(ieee754sp x) 38 { 39 return ieee754sp_class(x) >= IEEE754_CLASS_SNAN; 40 } 41 42 int ieee754sp_issnan(ieee754sp x) 43 { 44 assert(ieee754sp_isnan(x)); 45 return (SPMANT(x) & SP_MBIT(SP_MBITS-1)); 46 } 47 48 49 ieee754sp ieee754sp_xcpt(ieee754sp r, const char *op, ...) 50 { 51 struct ieee754xctx ax; 52 53 if (!TSTX()) 54 return r; 55 56 ax.op = op; 57 ax.rt = IEEE754_RT_SP; 58 ax.rv.sp = r; 59 va_start(ax.ap, op); 60 ieee754_xcpt(&ax); 61 return ax.rv.sp; 62 } 63 64 ieee754sp ieee754sp_nanxcpt(ieee754sp r, const char *op, ...) 65 { 66 struct ieee754xctx ax; 67 68 assert(ieee754sp_isnan(r)); 69 70 if (!ieee754sp_issnan(r)) /* QNAN does not cause invalid op !! */ 71 return r; 72 73 if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) { 74 /* not enabled convert to a quiet NaN */ 75 SPMANT(r) &= (~SP_MBIT(SP_MBITS-1)); 76 if (ieee754sp_isnan(r)) 77 return r; 78 else 79 return ieee754sp_indef(); 80 } 81 82 ax.op = op; 83 ax.rt = 0; 84 ax.rv.sp = r; 85 va_start(ax.ap, op); 86 ieee754_xcpt(&ax); 87 return ax.rv.sp; 88 } 89 90 ieee754sp ieee754sp_bestnan(ieee754sp x, ieee754sp y) 91 { 92 assert(ieee754sp_isnan(x)); 93 assert(ieee754sp_isnan(y)); 94 95 if (SPMANT(x) > SPMANT(y)) 96 return x; 97 else 98 return y; 99 } 100 101 102 static unsigned get_rounding(int sn, unsigned xm) 103 { 104 /* inexact must round of 3 bits 105 */ 106 if (xm & (SP_MBIT(3) - 1)) { 107 switch (ieee754_csr.rm) { 108 case IEEE754_RZ: 109 break; 110 case IEEE754_RN: 111 xm += 0x3 + ((xm >> 3) & 1); 112 /* xm += (xm&0x8)?0x4:0x3 */ 113 break; 114 case IEEE754_RU: /* toward +Infinity */ 115 if (!sn) /* ?? */ 116 xm += 0x8; 117 break; 118 case IEEE754_RD: /* toward -Infinity */ 119 if (sn) /* ?? */ 120 xm += 0x8; 121 break; 122 } 123 } 124 return xm; 125 } 126 127 128 /* generate a normal/denormal number with over,under handling 129 * sn is sign 130 * xe is an unbiased exponent 131 * xm is 3bit extended precision value. 132 */ 133 ieee754sp ieee754sp_format(int sn, int xe, unsigned xm) 134 { 135 assert(xm); /* we don't gen exact zeros (probably should) */ 136 137 assert((xm >> (SP_MBITS + 1 + 3)) == 0); /* no execess */ 138 assert(xm & (SP_HIDDEN_BIT << 3)); 139 140 if (xe < SP_EMIN) { 141 /* strip lower bits */ 142 int es = SP_EMIN - xe; 143 144 if (ieee754_csr.nod) { 145 SETCX(IEEE754_UNDERFLOW); 146 SETCX(IEEE754_INEXACT); 147 148 switch(ieee754_csr.rm) { 149 case IEEE754_RN: 150 return ieee754sp_zero(sn); 151 case IEEE754_RZ: 152 return ieee754sp_zero(sn); 153 case IEEE754_RU: /* toward +Infinity */ 154 if(sn == 0) 155 return ieee754sp_min(0); 156 else 157 return ieee754sp_zero(1); 158 case IEEE754_RD: /* toward -Infinity */ 159 if(sn == 0) 160 return ieee754sp_zero(0); 161 else 162 return ieee754sp_min(1); 163 } 164 } 165 166 if (xe == SP_EMIN - 1 167 && get_rounding(sn, xm) >> (SP_MBITS + 1 + 3)) 168 { 169 /* Not tiny after rounding */ 170 SETCX(IEEE754_INEXACT); 171 xm = get_rounding(sn, xm); 172 xm >>= 1; 173 /* Clear grs bits */ 174 xm &= ~(SP_MBIT(3) - 1); 175 xe++; 176 } 177 else { 178 /* sticky right shift es bits 179 */ 180 SPXSRSXn(es); 181 assert((xm & (SP_HIDDEN_BIT << 3)) == 0); 182 assert(xe == SP_EMIN); 183 } 184 } 185 if (xm & (SP_MBIT(3) - 1)) { 186 SETCX(IEEE754_INEXACT); 187 if ((xm & (SP_HIDDEN_BIT << 3)) == 0) { 188 SETCX(IEEE754_UNDERFLOW); 189 } 190 191 /* inexact must round of 3 bits 192 */ 193 xm = get_rounding(sn, xm); 194 /* adjust exponent for rounding add overflowing 195 */ 196 if (xm >> (SP_MBITS + 1 + 3)) { 197 /* add causes mantissa overflow */ 198 xm >>= 1; 199 xe++; 200 } 201 } 202 /* strip grs bits */ 203 xm >>= 3; 204 205 assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */ 206 assert(xe >= SP_EMIN); 207 208 if (xe > SP_EMAX) { 209 SETCX(IEEE754_OVERFLOW); 210 SETCX(IEEE754_INEXACT); 211 /* -O can be table indexed by (rm,sn) */ 212 switch (ieee754_csr.rm) { 213 case IEEE754_RN: 214 return ieee754sp_inf(sn); 215 case IEEE754_RZ: 216 return ieee754sp_max(sn); 217 case IEEE754_RU: /* toward +Infinity */ 218 if (sn == 0) 219 return ieee754sp_inf(0); 220 else 221 return ieee754sp_max(1); 222 case IEEE754_RD: /* toward -Infinity */ 223 if (sn == 0) 224 return ieee754sp_max(0); 225 else 226 return ieee754sp_inf(1); 227 } 228 } 229 /* gen norm/denorm/zero */ 230 231 if ((xm & SP_HIDDEN_BIT) == 0) { 232 /* we underflow (tiny/zero) */ 233 assert(xe == SP_EMIN); 234 if (ieee754_csr.mx & IEEE754_UNDERFLOW) 235 SETCX(IEEE754_UNDERFLOW); 236 return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm); 237 } else { 238 assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */ 239 assert(xm & SP_HIDDEN_BIT); 240 241 return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT); 242 } 243 } 244