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 * 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 "ieee754dp.h" 29 30 int ieee754dp_class(ieee754dp x) 31 { 32 COMPXDP; 33 EXPLODEXDP; 34 return xc; 35 } 36 37 int ieee754dp_isnan(ieee754dp x) 38 { 39 return ieee754dp_class(x) >= IEEE754_CLASS_SNAN; 40 } 41 42 int ieee754dp_issnan(ieee754dp x) 43 { 44 assert(ieee754dp_isnan(x)); 45 return ((DPMANT(x) & DP_MBIT(DP_MBITS-1)) == DP_MBIT(DP_MBITS-1)); 46 } 47 48 49 ieee754dp ieee754dp_xcpt(ieee754dp r, const char *op, ...) 50 { 51 struct ieee754xctx ax; 52 if (!TSTX()) 53 return r; 54 55 ax.op = op; 56 ax.rt = IEEE754_RT_DP; 57 ax.rv.dp = r; 58 va_start(ax.ap, op); 59 ieee754_xcpt(&ax); 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 return ax.rv.dp; 87 } 88 89 ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y) 90 { 91 assert(ieee754dp_isnan(x)); 92 assert(ieee754dp_isnan(y)); 93 94 if (DPMANT(x) > DPMANT(y)) 95 return x; 96 else 97 return y; 98 } 99 100 101 static u64 get_rounding(int sn, u64 xm) 102 { 103 /* inexact must round of 3 bits 104 */ 105 if (xm & (DP_MBIT(3) - 1)) { 106 switch (ieee754_csr.rm) { 107 case IEEE754_RZ: 108 break; 109 case IEEE754_RN: 110 xm += 0x3 + ((xm >> 3) & 1); 111 /* xm += (xm&0x8)?0x4:0x3 */ 112 break; 113 case IEEE754_RU: /* toward +Infinity */ 114 if (!sn) /* ?? */ 115 xm += 0x8; 116 break; 117 case IEEE754_RD: /* toward -Infinity */ 118 if (sn) /* ?? */ 119 xm += 0x8; 120 break; 121 } 122 } 123 return xm; 124 } 125 126 127 /* generate a normal/denormal number with over,under handling 128 * sn is sign 129 * xe is an unbiased exponent 130 * xm is 3bit extended precision value. 131 */ 132 ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 133 { 134 assert(xm); /* we don't gen exact zeros (probably should) */ 135 136 assert((xm >> (DP_MBITS + 1 + 3)) == 0); /* no execess */ 137 assert(xm & (DP_HIDDEN_BIT << 3)); 138 139 if (xe < DP_EMIN) { 140 /* strip lower bits */ 141 int es = DP_EMIN - xe; 142 143 if (ieee754_csr.nod) { 144 SETCX(IEEE754_UNDERFLOW); 145 SETCX(IEEE754_INEXACT); 146 147 switch(ieee754_csr.rm) { 148 case IEEE754_RN: 149 return ieee754dp_zero(sn); 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