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 va_end(ax.ap); 62 return ax.rv.sp; 63 } 64 65 ieee754sp ieee754sp_nanxcpt(ieee754sp r, const char *op, ...) 66 { 67 struct ieee754xctx ax; 68 69 assert(ieee754sp_isnan(r)); 70 71 if (!ieee754sp_issnan(r)) /* QNAN does not cause invalid op !! */ 72 return r; 73 74 if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) { 75 /* not enabled convert to a quiet NaN */ 76 SPMANT(r) &= (~SP_MBIT(SP_MBITS-1)); 77 if (ieee754sp_isnan(r)) 78 return r; 79 else 80 return ieee754sp_indef(); 81 } 82 83 ax.op = op; 84 ax.rt = 0; 85 ax.rv.sp = r; 86 va_start(ax.ap, op); 87 ieee754_xcpt(&ax); 88 va_end(ax.ap); 89 return ax.rv.sp; 90 } 91 92 ieee754sp ieee754sp_bestnan(ieee754sp x, ieee754sp y) 93 { 94 assert(ieee754sp_isnan(x)); 95 assert(ieee754sp_isnan(y)); 96 97 if (SPMANT(x) > SPMANT(y)) 98 return x; 99 else 100 return y; 101 } 102 103 104 static unsigned get_rounding(int sn, unsigned xm) 105 { 106 /* inexact must round of 3 bits 107 */ 108 if (xm & (SP_MBIT(3) - 1)) { 109 switch (ieee754_csr.rm) { 110 case IEEE754_RZ: 111 break; 112 case IEEE754_RN: 113 xm += 0x3 + ((xm >> 3) & 1); 114 /* xm += (xm&0x8)?0x4:0x3 */ 115 break; 116 case IEEE754_RU: /* toward +Infinity */ 117 if (!sn) /* ?? */ 118 xm += 0x8; 119 break; 120 case IEEE754_RD: /* toward -Infinity */ 121 if (sn) /* ?? */ 122 xm += 0x8; 123 break; 124 } 125 } 126 return xm; 127 } 128 129 130 /* generate a normal/denormal number with over,under handling 131 * sn is sign 132 * xe is an unbiased exponent 133 * xm is 3bit extended precision value. 134 */ 135 ieee754sp ieee754sp_format(int sn, int xe, unsigned xm) 136 { 137 assert(xm); /* we don't gen exact zeros (probably should) */ 138 139 assert((xm >> (SP_MBITS + 1 + 3)) == 0); /* no execess */ 140 assert(xm & (SP_HIDDEN_BIT << 3)); 141 142 if (xe < SP_EMIN) { 143 /* strip lower bits */ 144 int es = SP_EMIN - xe; 145 146 if (ieee754_csr.nod) { 147 SETCX(IEEE754_UNDERFLOW); 148 SETCX(IEEE754_INEXACT); 149 150 switch(ieee754_csr.rm) { 151 case IEEE754_RN: 152 return ieee754sp_zero(sn); 153 case IEEE754_RZ: 154 return ieee754sp_zero(sn); 155 case IEEE754_RU: /* toward +Infinity */ 156 if(sn == 0) 157 return ieee754sp_min(0); 158 else 159 return ieee754sp_zero(1); 160 case IEEE754_RD: /* toward -Infinity */ 161 if(sn == 0) 162 return ieee754sp_zero(0); 163 else 164 return ieee754sp_min(1); 165 } 166 } 167 168 if (xe == SP_EMIN - 1 169 && get_rounding(sn, xm) >> (SP_MBITS + 1 + 3)) 170 { 171 /* Not tiny after rounding */ 172 SETCX(IEEE754_INEXACT); 173 xm = get_rounding(sn, xm); 174 xm >>= 1; 175 /* Clear grs bits */ 176 xm &= ~(SP_MBIT(3) - 1); 177 xe++; 178 } 179 else { 180 /* sticky right shift es bits 181 */ 182 SPXSRSXn(es); 183 assert((xm & (SP_HIDDEN_BIT << 3)) == 0); 184 assert(xe == SP_EMIN); 185 } 186 } 187 if (xm & (SP_MBIT(3) - 1)) { 188 SETCX(IEEE754_INEXACT); 189 if ((xm & (SP_HIDDEN_BIT << 3)) == 0) { 190 SETCX(IEEE754_UNDERFLOW); 191 } 192 193 /* inexact must round of 3 bits 194 */ 195 xm = get_rounding(sn, xm); 196 /* adjust exponent for rounding add overflowing 197 */ 198 if (xm >> (SP_MBITS + 1 + 3)) { 199 /* add causes mantissa overflow */ 200 xm >>= 1; 201 xe++; 202 } 203 } 204 /* strip grs bits */ 205 xm >>= 3; 206 207 assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */ 208 assert(xe >= SP_EMIN); 209 210 if (xe > SP_EMAX) { 211 SETCX(IEEE754_OVERFLOW); 212 SETCX(IEEE754_INEXACT); 213 /* -O can be table indexed by (rm,sn) */ 214 switch (ieee754_csr.rm) { 215 case IEEE754_RN: 216 return ieee754sp_inf(sn); 217 case IEEE754_RZ: 218 return ieee754sp_max(sn); 219 case IEEE754_RU: /* toward +Infinity */ 220 if (sn == 0) 221 return ieee754sp_inf(0); 222 else 223 return ieee754sp_max(1); 224 case IEEE754_RD: /* toward -Infinity */ 225 if (sn == 0) 226 return ieee754sp_max(0); 227 else 228 return ieee754sp_inf(1); 229 } 230 } 231 /* gen norm/denorm/zero */ 232 233 if ((xm & SP_HIDDEN_BIT) == 0) { 234 /* we underflow (tiny/zero) */ 235 assert(xe == SP_EMIN); 236 if (ieee754_csr.mx & IEEE754_UNDERFLOW) 237 SETCX(IEEE754_UNDERFLOW); 238 return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm); 239 } else { 240 assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */ 241 assert(xm & SP_HIDDEN_BIT); 242 243 return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT); 244 } 245 } 246