1 /* IEEE754 floating point arithmetic 2 * double precision: common utilities 3 */ 4 /* 5 * MIPS floating point support 6 * Copyright (C) 1994-2000 Algorithmics Ltd. 7 * 8 * ######################################################################## 9 * 10 * This program is free software; you can distribute it and/or modify it 11 * under the terms of the GNU General Public License (Version 2) as 12 * published by the Free Software Foundation. 13 * 14 * This program is distributed in the hope it will be useful, but WITHOUT 15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 17 * for more details. 18 * 19 * You should have received a copy of the GNU General Public License along 20 * with this program; if not, write to the Free Software Foundation, Inc., 21 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. 22 * 23 * ######################################################################## 24 */ 25 26 27 #include "ieee754dp.h" 28 29 int ieee754dp_class(ieee754dp x) 30 { 31 COMPXDP; 32 EXPLODEXDP; 33 return xc; 34 } 35 36 int ieee754dp_isnan(ieee754dp x) 37 { 38 return ieee754dp_class(x) >= IEEE754_CLASS_SNAN; 39 } 40 41 int ieee754dp_issnan(ieee754dp x) 42 { 43 assert(ieee754dp_isnan(x)); 44 return ((DPMANT(x) & DP_MBIT(DP_MBITS-1)) == DP_MBIT(DP_MBITS-1)); 45 } 46 47 48 ieee754dp ieee754dp_xcpt(ieee754dp r, const char *op, ...) 49 { 50 struct ieee754xctx ax; 51 if (!TSTX()) 52 return r; 53 54 ax.op = op; 55 ax.rt = IEEE754_RT_DP; 56 ax.rv.dp = r; 57 va_start(ax.ap, op); 58 ieee754_xcpt(&ax); 59 va_end(ax.ap); 60 return ax.rv.dp; 61 } 62 63 ieee754dp ieee754dp_nanxcpt(ieee754dp r, const char *op, ...) 64 { 65 struct ieee754xctx ax; 66 67 assert(ieee754dp_isnan(r)); 68 69 if (!ieee754dp_issnan(r)) /* QNAN does not cause invalid op !! */ 70 return r; 71 72 if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) { 73 /* not enabled convert to a quiet NaN */ 74 DPMANT(r) &= (~DP_MBIT(DP_MBITS-1)); 75 if (ieee754dp_isnan(r)) 76 return r; 77 else 78 return ieee754dp_indef(); 79 } 80 81 ax.op = op; 82 ax.rt = 0; 83 ax.rv.dp = r; 84 va_start(ax.ap, op); 85 ieee754_xcpt(&ax); 86 va_end(ax.ap); 87 return ax.rv.dp; 88 } 89 90 ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y) 91 { 92 assert(ieee754dp_isnan(x)); 93 assert(ieee754dp_isnan(y)); 94 95 if (DPMANT(x) > DPMANT(y)) 96 return x; 97 else 98 return y; 99 } 100 101 102 static u64 get_rounding(int sn, u64 xm) 103 { 104 /* inexact must round of 3 bits 105 */ 106 if (xm & (DP_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 ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 134 { 135 assert(xm); /* we don't gen exact zeros (probably should) */ 136 137 assert((xm >> (DP_MBITS + 1 + 3)) == 0); /* no execess */ 138 assert(xm & (DP_HIDDEN_BIT << 3)); 139 140 if (xe < DP_EMIN) { 141 /* strip lower bits */ 142 int es = DP_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 case IEEE754_RZ: 151 return ieee754dp_zero(sn); 152 case IEEE754_RU: /* toward +Infinity */ 153 if(sn == 0) 154 return ieee754dp_min(0); 155 else 156 return ieee754dp_zero(1); 157 case IEEE754_RD: /* toward -Infinity */ 158 if(sn == 0) 159 return ieee754dp_zero(0); 160 else 161 return ieee754dp_min(1); 162 } 163 } 164 165 if (xe == DP_EMIN - 1 166 && get_rounding(sn, xm) >> (DP_MBITS + 1 + 3)) 167 { 168 /* Not tiny after rounding */ 169 SETCX(IEEE754_INEXACT); 170 xm = get_rounding(sn, xm); 171 xm >>= 1; 172 /* Clear grs bits */ 173 xm &= ~(DP_MBIT(3) - 1); 174 xe++; 175 } 176 else { 177 /* sticky right shift es bits 178 */ 179 xm = XDPSRS(xm, es); 180 xe += es; 181 assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 182 assert(xe == DP_EMIN); 183 } 184 } 185 if (xm & (DP_MBIT(3) - 1)) { 186 SETCX(IEEE754_INEXACT); 187 if ((xm & (DP_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 >> (DP_MBITS + 3 + 1)) { 197 /* add causes mantissa overflow */ 198 xm >>= 1; 199 xe++; 200 } 201 } 202 /* strip grs bits */ 203 xm >>= 3; 204 205 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ 206 assert(xe >= DP_EMIN); 207 208 if (xe > DP_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 ieee754dp_inf(sn); 215 case IEEE754_RZ: 216 return ieee754dp_max(sn); 217 case IEEE754_RU: /* toward +Infinity */ 218 if (sn == 0) 219 return ieee754dp_inf(0); 220 else 221 return ieee754dp_max(1); 222 case IEEE754_RD: /* toward -Infinity */ 223 if (sn == 0) 224 return ieee754dp_max(0); 225 else 226 return ieee754dp_inf(1); 227 } 228 } 229 /* gen norm/denorm/zero */ 230 231 if ((xm & DP_HIDDEN_BIT) == 0) { 232 /* we underflow (tiny/zero) */ 233 assert(xe == DP_EMIN); 234 if (ieee754_csr.mx & IEEE754_UNDERFLOW) 235 SETCX(IEEE754_UNDERFLOW); 236 return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 237 } else { 238 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ 239 assert(xm & DP_HIDDEN_BIT); 240 241 return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 242 } 243 } 244