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 static inline int ieee754dp_isnan(union ieee754dp x) 34 { 35 return ieee754_class_nan(ieee754dp_class(x)); 36 } 37 38 static inline int ieee754dp_issnan(union ieee754dp x) 39 { 40 int qbit; 41 42 assert(ieee754dp_isnan(x)); 43 qbit = (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_FBITS - 1); 44 return ieee754_csr.nan2008 ^ qbit; 45 } 46 47 48 /* 49 * Raise the Invalid Operation IEEE 754 exception 50 * and convert the signaling NaN supplied to a quiet NaN. 51 */ 52 union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r) 53 { 54 assert(ieee754dp_issnan(r)); 55 56 ieee754_setcx(IEEE754_INVALID_OPERATION); 57 if (ieee754_csr.nan2008) { 58 DPMANT(r) |= DP_MBIT(DP_FBITS - 1); 59 } else { 60 DPMANT(r) &= ~DP_MBIT(DP_FBITS - 1); 61 if (!ieee754dp_isnan(r)) 62 DPMANT(r) |= DP_MBIT(DP_FBITS - 2); 63 } 64 65 return r; 66 } 67 68 static u64 ieee754dp_get_rounding(int sn, u64 xm) 69 { 70 /* inexact must round of 3 bits 71 */ 72 if (xm & (DP_MBIT(3) - 1)) { 73 switch (ieee754_csr.rm) { 74 case FPU_CSR_RZ: 75 break; 76 case FPU_CSR_RN: 77 xm += 0x3 + ((xm >> 3) & 1); 78 /* xm += (xm&0x8)?0x4:0x3 */ 79 break; 80 case FPU_CSR_RU: /* toward +Infinity */ 81 if (!sn) /* ?? */ 82 xm += 0x8; 83 break; 84 case FPU_CSR_RD: /* toward -Infinity */ 85 if (sn) /* ?? */ 86 xm += 0x8; 87 break; 88 } 89 } 90 return xm; 91 } 92 93 94 /* generate a normal/denormal number with over,under handling 95 * sn is sign 96 * xe is an unbiased exponent 97 * xm is 3bit extended precision value. 98 */ 99 union ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 100 { 101 assert(xm); /* we don't gen exact zeros (probably should) */ 102 103 assert((xm >> (DP_FBITS + 1 + 3)) == 0); /* no excess */ 104 assert(xm & (DP_HIDDEN_BIT << 3)); 105 106 if (xe < DP_EMIN) { 107 /* strip lower bits */ 108 int es = DP_EMIN - xe; 109 110 if (ieee754_csr.nod) { 111 ieee754_setcx(IEEE754_UNDERFLOW); 112 ieee754_setcx(IEEE754_INEXACT); 113 114 switch(ieee754_csr.rm) { 115 case FPU_CSR_RN: 116 case FPU_CSR_RZ: 117 return ieee754dp_zero(sn); 118 case FPU_CSR_RU: /* toward +Infinity */ 119 if (sn == 0) 120 return ieee754dp_min(0); 121 else 122 return ieee754dp_zero(1); 123 case FPU_CSR_RD: /* toward -Infinity */ 124 if (sn == 0) 125 return ieee754dp_zero(0); 126 else 127 return ieee754dp_min(1); 128 } 129 } 130 131 if (xe == DP_EMIN - 1 && 132 ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3)) 133 { 134 /* Not tiny after rounding */ 135 ieee754_setcx(IEEE754_INEXACT); 136 xm = ieee754dp_get_rounding(sn, xm); 137 xm >>= 1; 138 /* Clear grs bits */ 139 xm &= ~(DP_MBIT(3) - 1); 140 xe++; 141 } 142 else { 143 /* sticky right shift es bits 144 */ 145 xm = XDPSRS(xm, es); 146 xe += es; 147 assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 148 assert(xe == DP_EMIN); 149 } 150 } 151 if (xm & (DP_MBIT(3) - 1)) { 152 ieee754_setcx(IEEE754_INEXACT); 153 if ((xm & (DP_HIDDEN_BIT << 3)) == 0) { 154 ieee754_setcx(IEEE754_UNDERFLOW); 155 } 156 157 /* inexact must round of 3 bits 158 */ 159 xm = ieee754dp_get_rounding(sn, xm); 160 /* adjust exponent for rounding add overflowing 161 */ 162 if (xm >> (DP_FBITS + 3 + 1)) { 163 /* add causes mantissa overflow */ 164 xm >>= 1; 165 xe++; 166 } 167 } 168 /* strip grs bits */ 169 xm >>= 3; 170 171 assert((xm >> (DP_FBITS + 1)) == 0); /* no excess */ 172 assert(xe >= DP_EMIN); 173 174 if (xe > DP_EMAX) { 175 ieee754_setcx(IEEE754_OVERFLOW); 176 ieee754_setcx(IEEE754_INEXACT); 177 /* -O can be table indexed by (rm,sn) */ 178 switch (ieee754_csr.rm) { 179 case FPU_CSR_RN: 180 return ieee754dp_inf(sn); 181 case FPU_CSR_RZ: 182 return ieee754dp_max(sn); 183 case FPU_CSR_RU: /* toward +Infinity */ 184 if (sn == 0) 185 return ieee754dp_inf(0); 186 else 187 return ieee754dp_max(1); 188 case FPU_CSR_RD: /* toward -Infinity */ 189 if (sn == 0) 190 return ieee754dp_max(0); 191 else 192 return ieee754dp_inf(1); 193 } 194 } 195 /* gen norm/denorm/zero */ 196 197 if ((xm & DP_HIDDEN_BIT) == 0) { 198 /* we underflow (tiny/zero) */ 199 assert(xe == DP_EMIN); 200 if (ieee754_csr.mx & IEEE754_UNDERFLOW) 201 ieee754_setcx(IEEE754_UNDERFLOW); 202 return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 203 } else { 204 assert((xm >> (DP_FBITS + 1)) == 0); /* no excess */ 205 assert(xm & DP_HIDDEN_BIT); 206 207 return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 208 } 209 } 210