1 /* 2 * Ported from a work by Andreas Grabher for Previous, NeXT Computer Emulator, 3 * derived from NetBSD M68040 FPSP functions, 4 * derived from release 2a of the SoftFloat IEC/IEEE Floating-point Arithmetic 5 * Package. Those parts of the code (and some later contributions) are 6 * provided under that license, as detailed below. 7 * It has subsequently been modified by contributors to the QEMU Project, 8 * so some portions are provided under: 9 * the SoftFloat-2a license 10 * the BSD license 11 * GPL-v2-or-later 12 * 13 * Any future contributions to this file will be taken to be licensed under 14 * the Softfloat-2a license unless specifically indicated otherwise. 15 */ 16 17 /* Portions of this work are licensed under the terms of the GNU GPL, 18 * version 2 or later. See the COPYING file in the top-level directory. 19 */ 20 21 #include "qemu/osdep.h" 22 #include "softfloat.h" 23 #include "fpu/softfloat-macros.h" 24 25 static floatx80 propagateFloatx80NaNOneArg(floatx80 a, float_status *status) 26 { 27 if (floatx80_is_signaling_nan(a, status)) { 28 float_raise(float_flag_invalid, status); 29 } 30 31 if (status->default_nan_mode) { 32 return floatx80_default_nan(status); 33 } 34 35 return floatx80_maybe_silence_nan(a, status); 36 } 37 38 /*---------------------------------------------------------------------------- 39 | Returns the modulo remainder of the extended double-precision floating-point 40 | value `a' with respect to the corresponding value `b'. 41 *----------------------------------------------------------------------------*/ 42 43 floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status) 44 { 45 flag aSign, zSign; 46 int32_t aExp, bExp, expDiff; 47 uint64_t aSig0, aSig1, bSig; 48 uint64_t qTemp, term0, term1; 49 50 aSig0 = extractFloatx80Frac(a); 51 aExp = extractFloatx80Exp(a); 52 aSign = extractFloatx80Sign(a); 53 bSig = extractFloatx80Frac(b); 54 bExp = extractFloatx80Exp(b); 55 56 if (aExp == 0x7FFF) { 57 if ((uint64_t) (aSig0 << 1) 58 || ((bExp == 0x7FFF) && (uint64_t) (bSig << 1))) { 59 return propagateFloatx80NaN(a, b, status); 60 } 61 goto invalid; 62 } 63 if (bExp == 0x7FFF) { 64 if ((uint64_t) (bSig << 1)) { 65 return propagateFloatx80NaN(a, b, status); 66 } 67 return a; 68 } 69 if (bExp == 0) { 70 if (bSig == 0) { 71 invalid: 72 float_raise(float_flag_invalid, status); 73 return floatx80_default_nan(status); 74 } 75 normalizeFloatx80Subnormal(bSig, &bExp, &bSig); 76 } 77 if (aExp == 0) { 78 if ((uint64_t) (aSig0 << 1) == 0) { 79 return a; 80 } 81 normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0); 82 } 83 bSig |= LIT64(0x8000000000000000); 84 zSign = aSign; 85 expDiff = aExp - bExp; 86 aSig1 = 0; 87 if (expDiff < 0) { 88 return a; 89 } 90 qTemp = (bSig <= aSig0); 91 if (qTemp) { 92 aSig0 -= bSig; 93 } 94 expDiff -= 64; 95 while (0 < expDiff) { 96 qTemp = estimateDiv128To64(aSig0, aSig1, bSig); 97 qTemp = (2 < qTemp) ? qTemp - 2 : 0; 98 mul64To128(bSig, qTemp, &term0, &term1); 99 sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 100 shortShift128Left(aSig0, aSig1, 62, &aSig0, &aSig1); 101 } 102 expDiff += 64; 103 if (0 < expDiff) { 104 qTemp = estimateDiv128To64(aSig0, aSig1, bSig); 105 qTemp = (2 < qTemp) ? qTemp - 2 : 0; 106 qTemp >>= 64 - expDiff; 107 mul64To128(bSig, qTemp << (64 - expDiff), &term0, &term1); 108 sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 109 shortShift128Left(0, bSig, 64 - expDiff, &term0, &term1); 110 while (le128(term0, term1, aSig0, aSig1)) { 111 ++qTemp; 112 sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 113 } 114 } 115 return 116 normalizeRoundAndPackFloatx80( 117 80, zSign, bExp + expDiff, aSig0, aSig1, status); 118 } 119 120 /*---------------------------------------------------------------------------- 121 | Returns the mantissa of the extended double-precision floating-point 122 | value `a'. 123 *----------------------------------------------------------------------------*/ 124 125 floatx80 floatx80_getman(floatx80 a, float_status *status) 126 { 127 flag aSign; 128 int32_t aExp; 129 uint64_t aSig; 130 131 aSig = extractFloatx80Frac(a); 132 aExp = extractFloatx80Exp(a); 133 aSign = extractFloatx80Sign(a); 134 135 if (aExp == 0x7FFF) { 136 if ((uint64_t) (aSig << 1)) { 137 return propagateFloatx80NaNOneArg(a , status); 138 } 139 float_raise(float_flag_invalid , status); 140 return floatx80_default_nan(status); 141 } 142 143 if (aExp == 0) { 144 if (aSig == 0) { 145 return packFloatx80(aSign, 0, 0); 146 } 147 normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 148 } 149 150 return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign, 151 0x3FFF, aSig, 0, status); 152 } 153 154 /*---------------------------------------------------------------------------- 155 | Returns the exponent of the extended double-precision floating-point 156 | value `a' as an extended double-precision value. 157 *----------------------------------------------------------------------------*/ 158 159 floatx80 floatx80_getexp(floatx80 a, float_status *status) 160 { 161 flag aSign; 162 int32_t aExp; 163 uint64_t aSig; 164 165 aSig = extractFloatx80Frac(a); 166 aExp = extractFloatx80Exp(a); 167 aSign = extractFloatx80Sign(a); 168 169 if (aExp == 0x7FFF) { 170 if ((uint64_t) (aSig << 1)) { 171 return propagateFloatx80NaNOneArg(a , status); 172 } 173 float_raise(float_flag_invalid , status); 174 return floatx80_default_nan(status); 175 } 176 177 if (aExp == 0) { 178 if (aSig == 0) { 179 return packFloatx80(aSign, 0, 0); 180 } 181 normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 182 } 183 184 return int32_to_floatx80(aExp - 0x3FFF, status); 185 } 186 187 /*---------------------------------------------------------------------------- 188 | Scales extended double-precision floating-point value in operand `a' by 189 | value `b'. The function truncates the value in the second operand 'b' to 190 | an integral value and adds that value to the exponent of the operand 'a'. 191 | The operation performed according to the IEC/IEEE Standard for Binary 192 | Floating-Point Arithmetic. 193 *----------------------------------------------------------------------------*/ 194 195 floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status) 196 { 197 flag aSign, bSign; 198 int32_t aExp, bExp, shiftCount; 199 uint64_t aSig, bSig; 200 201 aSig = extractFloatx80Frac(a); 202 aExp = extractFloatx80Exp(a); 203 aSign = extractFloatx80Sign(a); 204 bSig = extractFloatx80Frac(b); 205 bExp = extractFloatx80Exp(b); 206 bSign = extractFloatx80Sign(b); 207 208 if (bExp == 0x7FFF) { 209 if ((uint64_t) (bSig << 1) || 210 ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) { 211 return propagateFloatx80NaN(a, b, status); 212 } 213 float_raise(float_flag_invalid , status); 214 return floatx80_default_nan(status); 215 } 216 if (aExp == 0x7FFF) { 217 if ((uint64_t) (aSig << 1)) { 218 return propagateFloatx80NaN(a, b, status); 219 } 220 return packFloatx80(aSign, floatx80_infinity.high, 221 floatx80_infinity.low); 222 } 223 if (aExp == 0) { 224 if (aSig == 0) { 225 return packFloatx80(aSign, 0, 0); 226 } 227 if (bExp < 0x3FFF) { 228 return a; 229 } 230 normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 231 } 232 233 if (bExp < 0x3FFF) { 234 return a; 235 } 236 237 if (0x400F < bExp) { 238 aExp = bSign ? -0x6001 : 0xE000; 239 return roundAndPackFloatx80(status->floatx80_rounding_precision, 240 aSign, aExp, aSig, 0, status); 241 } 242 243 shiftCount = 0x403E - bExp; 244 bSig >>= shiftCount; 245 aExp = bSign ? (aExp - bSig) : (aExp + bSig); 246 247 return roundAndPackFloatx80(status->floatx80_rounding_precision, 248 aSign, aExp, aSig, 0, status); 249 } 250