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 * This program is free software; you can distribute it and/or modify it 9 * under the terms of the GNU General Public License (Version 2) as 10 * published by the Free Software Foundation. 11 * 12 * This program is distributed in the hope it will be useful, but WITHOUT 13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 14 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 15 * for more details. 16 * 17 * You should have received a copy of the GNU General Public License along 18 * with this program; if not, write to the Free Software Foundation, Inc., 19 * 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 20 */ 21 22 #include <linux/compiler.h> 23 24 #include "ieee754dp.h" 25 26 int ieee754dp_class(union ieee754dp x) 27 { 28 COMPXDP; 29 EXPLODEXDP; 30 return xc; 31 } 32 33 int ieee754dp_isnan(union ieee754dp x) 34 { 35 return ieee754dp_class(x) >= IEEE754_CLASS_SNAN; 36 } 37 38 static inline int ieee754dp_issnan(union ieee754dp x) 39 { 40 assert(ieee754dp_isnan(x)); 41 return (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_FBITS - 1); 42 } 43 44 45 union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r) 46 { 47 assert(ieee754dp_isnan(r)); 48 49 if (!ieee754dp_issnan(r)) /* QNAN does not cause invalid op !! */ 50 return r; 51 52 if (!ieee754_setandtestcx(IEEE754_INVALID_OPERATION)) { 53 /* not enabled convert to a quiet NaN */ 54 DPMANT(r) &= (~DP_MBIT(DP_FBITS-1)); 55 if (ieee754dp_isnan(r)) 56 return r; 57 else 58 return ieee754dp_indef(); 59 } 60 61 return r; 62 } 63 64 static u64 ieee754dp_get_rounding(int sn, u64 xm) 65 { 66 /* inexact must round of 3 bits 67 */ 68 if (xm & (DP_MBIT(3) - 1)) { 69 switch (ieee754_csr.rm) { 70 case FPU_CSR_RZ: 71 break; 72 case FPU_CSR_RN: 73 xm += 0x3 + ((xm >> 3) & 1); 74 /* xm += (xm&0x8)?0x4:0x3 */ 75 break; 76 case FPU_CSR_RU: /* toward +Infinity */ 77 if (!sn) /* ?? */ 78 xm += 0x8; 79 break; 80 case FPU_CSR_RD: /* toward -Infinity */ 81 if (sn) /* ?? */ 82 xm += 0x8; 83 break; 84 } 85 } 86 return xm; 87 } 88 89 90 /* generate a normal/denormal number with over,under handling 91 * sn is sign 92 * xe is an unbiased exponent 93 * xm is 3bit extended precision value. 94 */ 95 union ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 96 { 97 assert(xm); /* we don't gen exact zeros (probably should) */ 98 99 assert((xm >> (DP_FBITS + 1 + 3)) == 0); /* no execess */ 100 assert(xm & (DP_HIDDEN_BIT << 3)); 101 102 if (xe < DP_EMIN) { 103 /* strip lower bits */ 104 int es = DP_EMIN - xe; 105 106 if (ieee754_csr.nod) { 107 ieee754_setcx(IEEE754_UNDERFLOW); 108 ieee754_setcx(IEEE754_INEXACT); 109 110 switch(ieee754_csr.rm) { 111 case FPU_CSR_RN: 112 case FPU_CSR_RZ: 113 return ieee754dp_zero(sn); 114 case FPU_CSR_RU: /* toward +Infinity */ 115 if (sn == 0) 116 return ieee754dp_min(0); 117 else 118 return ieee754dp_zero(1); 119 case FPU_CSR_RD: /* toward -Infinity */ 120 if (sn == 0) 121 return ieee754dp_zero(0); 122 else 123 return ieee754dp_min(1); 124 } 125 } 126 127 if (xe == DP_EMIN - 1 && 128 ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3)) 129 { 130 /* Not tiny after rounding */ 131 ieee754_setcx(IEEE754_INEXACT); 132 xm = ieee754dp_get_rounding(sn, xm); 133 xm >>= 1; 134 /* Clear grs bits */ 135 xm &= ~(DP_MBIT(3) - 1); 136 xe++; 137 } 138 else { 139 /* sticky right shift es bits 140 */ 141 xm = XDPSRS(xm, es); 142 xe += es; 143 assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 144 assert(xe == DP_EMIN); 145 } 146 } 147 if (xm & (DP_MBIT(3) - 1)) { 148 ieee754_setcx(IEEE754_INEXACT); 149 if ((xm & (DP_HIDDEN_BIT << 3)) == 0) { 150 ieee754_setcx(IEEE754_UNDERFLOW); 151 } 152 153 /* inexact must round of 3 bits 154 */ 155 xm = ieee754dp_get_rounding(sn, xm); 156 /* adjust exponent for rounding add overflowing 157 */ 158 if (xm >> (DP_FBITS + 3 + 1)) { 159 /* add causes mantissa overflow */ 160 xm >>= 1; 161 xe++; 162 } 163 } 164 /* strip grs bits */ 165 xm >>= 3; 166 167 assert((xm >> (DP_FBITS + 1)) == 0); /* no execess */ 168 assert(xe >= DP_EMIN); 169 170 if (xe > DP_EMAX) { 171 ieee754_setcx(IEEE754_OVERFLOW); 172 ieee754_setcx(IEEE754_INEXACT); 173 /* -O can be table indexed by (rm,sn) */ 174 switch (ieee754_csr.rm) { 175 case FPU_CSR_RN: 176 return ieee754dp_inf(sn); 177 case FPU_CSR_RZ: 178 return ieee754dp_max(sn); 179 case FPU_CSR_RU: /* toward +Infinity */ 180 if (sn == 0) 181 return ieee754dp_inf(0); 182 else 183 return ieee754dp_max(1); 184 case FPU_CSR_RD: /* toward -Infinity */ 185 if (sn == 0) 186 return ieee754dp_max(0); 187 else 188 return ieee754dp_inf(1); 189 } 190 } 191 /* gen norm/denorm/zero */ 192 193 if ((xm & DP_HIDDEN_BIT) == 0) { 194 /* we underflow (tiny/zero) */ 195 assert(xe == DP_EMIN); 196 if (ieee754_csr.mx & IEEE754_UNDERFLOW) 197 ieee754_setcx(IEEE754_UNDERFLOW); 198 return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 199 } else { 200 assert((xm >> (DP_FBITS + 1)) == 0); /* no execess */ 201 assert(xm & DP_HIDDEN_BIT); 202 203 return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 204 } 205 } 206