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 return ieee754dp_zero(sn); 152 case IEEE754_RZ: 153 return ieee754dp_zero(sn); 154 case IEEE754_RU: /* toward +Infinity */ 155 if(sn == 0) 156 return ieee754dp_min(0); 157 else 158 return ieee754dp_zero(1); 159 case IEEE754_RD: /* toward -Infinity */ 160 if(sn == 0) 161 return ieee754dp_zero(0); 162 else 163 return ieee754dp_min(1); 164 } 165 } 166 167 if (xe == DP_EMIN - 1 168 && get_rounding(sn, xm) >> (DP_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 &= ~(DP_MBIT(3) - 1); 176 xe++; 177 } 178 else { 179 /* sticky right shift es bits 180 */ 181 xm = XDPSRS(xm, es); 182 xe += es; 183 assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 184 assert(xe == DP_EMIN); 185 } 186 } 187 if (xm & (DP_MBIT(3) - 1)) { 188 SETCX(IEEE754_INEXACT); 189 if ((xm & (DP_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 >> (DP_MBITS + 3 + 1)) { 199 /* add causes mantissa overflow */ 200 xm >>= 1; 201 xe++; 202 } 203 } 204 /* strip grs bits */ 205 xm >>= 3; 206 207 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ 208 assert(xe >= DP_EMIN); 209 210 if (xe > DP_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 ieee754dp_inf(sn); 217 case IEEE754_RZ: 218 return ieee754dp_max(sn); 219 case IEEE754_RU: /* toward +Infinity */ 220 if (sn == 0) 221 return ieee754dp_inf(0); 222 else 223 return ieee754dp_max(1); 224 case IEEE754_RD: /* toward -Infinity */ 225 if (sn == 0) 226 return ieee754dp_max(0); 227 else 228 return ieee754dp_inf(1); 229 } 230 } 231 /* gen norm/denorm/zero */ 232 233 if ((xm & DP_HIDDEN_BIT) == 0) { 234 /* we underflow (tiny/zero) */ 235 assert(xe == DP_EMIN); 236 if (ieee754_csr.mx & IEEE754_UNDERFLOW) 237 SETCX(IEEE754_UNDERFLOW); 238 return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 239 } else { 240 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ 241 assert(xm & DP_HIDDEN_BIT); 242 243 return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 244 } 245 } 246