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 va_end(ax.ap); 61 return ax.rv.dp; 62 } 63 64 ieee754dp ieee754dp_nanxcpt(ieee754dp r, const char *op, ...) 65 { 66 struct ieee754xctx ax; 67 68 assert(ieee754dp_isnan(r)); 69 70 if (!ieee754dp_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 DPMANT(r) &= (~DP_MBIT(DP_MBITS-1)); 76 if (ieee754dp_isnan(r)) 77 return r; 78 else 79 return ieee754dp_indef(); 80 } 81 82 ax.op = op; 83 ax.rt = 0; 84 ax.rv.dp = r; 85 va_start(ax.ap, op); 86 ieee754_xcpt(&ax); 87 va_end(ax.ap); 88 return ax.rv.dp; 89 } 90 91 ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y) 92 { 93 assert(ieee754dp_isnan(x)); 94 assert(ieee754dp_isnan(y)); 95 96 if (DPMANT(x) > DPMANT(y)) 97 return x; 98 else 99 return y; 100 } 101 102 103 static u64 get_rounding(int sn, u64 xm) 104 { 105 /* inexact must round of 3 bits 106 */ 107 if (xm & (DP_MBIT(3) - 1)) { 108 switch (ieee754_csr.rm) { 109 case IEEE754_RZ: 110 break; 111 case IEEE754_RN: 112 xm += 0x3 + ((xm >> 3) & 1); 113 /* xm += (xm&0x8)?0x4:0x3 */ 114 break; 115 case IEEE754_RU: /* toward +Infinity */ 116 if (!sn) /* ?? */ 117 xm += 0x8; 118 break; 119 case IEEE754_RD: /* toward -Infinity */ 120 if (sn) /* ?? */ 121 xm += 0x8; 122 break; 123 } 124 } 125 return xm; 126 } 127 128 129 /* generate a normal/denormal number with over,under handling 130 * sn is sign 131 * xe is an unbiased exponent 132 * xm is 3bit extended precision value. 133 */ 134 ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 135 { 136 assert(xm); /* we don't gen exact zeros (probably should) */ 137 138 assert((xm >> (DP_MBITS + 1 + 3)) == 0); /* no execess */ 139 assert(xm & (DP_HIDDEN_BIT << 3)); 140 141 if (xe < DP_EMIN) { 142 /* strip lower bits */ 143 int es = DP_EMIN - xe; 144 145 if (ieee754_csr.nod) { 146 SETCX(IEEE754_UNDERFLOW); 147 SETCX(IEEE754_INEXACT); 148 149 switch(ieee754_csr.rm) { 150 case IEEE754_RN: 151 case IEEE754_RZ: 152 return ieee754dp_zero(sn); 153 case IEEE754_RU: /* toward +Infinity */ 154 if(sn == 0) 155 return ieee754dp_min(0); 156 else 157 return ieee754dp_zero(1); 158 case IEEE754_RD: /* toward -Infinity */ 159 if(sn == 0) 160 return ieee754dp_zero(0); 161 else 162 return ieee754dp_min(1); 163 } 164 } 165 166 if (xe == DP_EMIN - 1 167 && get_rounding(sn, xm) >> (DP_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 &= ~(DP_MBIT(3) - 1); 175 xe++; 176 } 177 else { 178 /* sticky right shift es bits 179 */ 180 xm = XDPSRS(xm, es); 181 xe += es; 182 assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 183 assert(xe == DP_EMIN); 184 } 185 } 186 if (xm & (DP_MBIT(3) - 1)) { 187 SETCX(IEEE754_INEXACT); 188 if ((xm & (DP_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 >> (DP_MBITS + 3 + 1)) { 198 /* add causes mantissa overflow */ 199 xm >>= 1; 200 xe++; 201 } 202 } 203 /* strip grs bits */ 204 xm >>= 3; 205 206 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ 207 assert(xe >= DP_EMIN); 208 209 if (xe > DP_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 ieee754dp_inf(sn); 216 case IEEE754_RZ: 217 return ieee754dp_max(sn); 218 case IEEE754_RU: /* toward +Infinity */ 219 if (sn == 0) 220 return ieee754dp_inf(0); 221 else 222 return ieee754dp_max(1); 223 case IEEE754_RD: /* toward -Infinity */ 224 if (sn == 0) 225 return ieee754dp_max(0); 226 else 227 return ieee754dp_inf(1); 228 } 229 } 230 /* gen norm/denorm/zero */ 231 232 if ((xm & DP_HIDDEN_BIT) == 0) { 233 /* we underflow (tiny/zero) */ 234 assert(xe == DP_EMIN); 235 if (ieee754_csr.mx & IEEE754_UNDERFLOW) 236 SETCX(IEEE754_UNDERFLOW); 237 return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 238 } else { 239 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ 240 assert(xm & DP_HIDDEN_BIT); 241 242 return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 243 } 244 } 245