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 case IEEE754_RZ: 153 return ieee754sp_zero(sn); 154 case IEEE754_RU: /* toward +Infinity */ 155 if(sn == 0) 156 return ieee754sp_min(0); 157 else 158 return ieee754sp_zero(1); 159 case IEEE754_RD: /* toward -Infinity */ 160 if(sn == 0) 161 return ieee754sp_zero(0); 162 else 163 return ieee754sp_min(1); 164 } 165 } 166 167 if (xe == SP_EMIN - 1 168 && get_rounding(sn, xm) >> (SP_MBITS + 1 + 3)) 169 { 170 /* Not tiny after rounding */ 171 SETCX(IEEE754_INEXACT); 172 xm = get_rounding(sn, xm); 173 xm >>= 1; 174 /* Clear grs bits */ 175 xm &= ~(SP_MBIT(3) - 1); 176 xe++; 177 } 178 else { 179 /* sticky right shift es bits 180 */ 181 SPXSRSXn(es); 182 assert((xm & (SP_HIDDEN_BIT << 3)) == 0); 183 assert(xe == SP_EMIN); 184 } 185 } 186 if (xm & (SP_MBIT(3) - 1)) { 187 SETCX(IEEE754_INEXACT); 188 if ((xm & (SP_HIDDEN_BIT << 3)) == 0) { 189 SETCX(IEEE754_UNDERFLOW); 190 } 191 192 /* inexact must round of 3 bits 193 */ 194 xm = get_rounding(sn, xm); 195 /* adjust exponent for rounding add overflowing 196 */ 197 if (xm >> (SP_MBITS + 1 + 3)) { 198 /* add causes mantissa overflow */ 199 xm >>= 1; 200 xe++; 201 } 202 } 203 /* strip grs bits */ 204 xm >>= 3; 205 206 assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */ 207 assert(xe >= SP_EMIN); 208 209 if (xe > SP_EMAX) { 210 SETCX(IEEE754_OVERFLOW); 211 SETCX(IEEE754_INEXACT); 212 /* -O can be table indexed by (rm,sn) */ 213 switch (ieee754_csr.rm) { 214 case IEEE754_RN: 215 return ieee754sp_inf(sn); 216 case IEEE754_RZ: 217 return ieee754sp_max(sn); 218 case IEEE754_RU: /* toward +Infinity */ 219 if (sn == 0) 220 return ieee754sp_inf(0); 221 else 222 return ieee754sp_max(1); 223 case IEEE754_RD: /* toward -Infinity */ 224 if (sn == 0) 225 return ieee754sp_max(0); 226 else 227 return ieee754sp_inf(1); 228 } 229 } 230 /* gen norm/denorm/zero */ 231 232 if ((xm & SP_HIDDEN_BIT) == 0) { 233 /* we underflow (tiny/zero) */ 234 assert(xe == SP_EMIN); 235 if (ieee754_csr.mx & IEEE754_UNDERFLOW) 236 SETCX(IEEE754_UNDERFLOW); 237 return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm); 238 } else { 239 assert((xm >> (SP_MBITS + 1)) == 0); /* no execess */ 240 assert(xm & SP_HIDDEN_BIT); 241 242 return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT); 243 } 244 } 245