1591596b7SLaurent Vivier /* 2591596b7SLaurent Vivier * Ported from a work by Andreas Grabher for Previous, NeXT Computer Emulator, 3591596b7SLaurent Vivier * derived from NetBSD M68040 FPSP functions, 4591596b7SLaurent Vivier * derived from release 2a of the SoftFloat IEC/IEEE Floating-point Arithmetic 5591596b7SLaurent Vivier * Package. Those parts of the code (and some later contributions) are 6591596b7SLaurent Vivier * provided under that license, as detailed below. 7591596b7SLaurent Vivier * It has subsequently been modified by contributors to the QEMU Project, 8591596b7SLaurent Vivier * so some portions are provided under: 9591596b7SLaurent Vivier * the SoftFloat-2a license 10591596b7SLaurent Vivier * the BSD license 11591596b7SLaurent Vivier * GPL-v2-or-later 12591596b7SLaurent Vivier * 13591596b7SLaurent Vivier * Any future contributions to this file will be taken to be licensed under 14591596b7SLaurent Vivier * the Softfloat-2a license unless specifically indicated otherwise. 15591596b7SLaurent Vivier */ 16591596b7SLaurent Vivier 17591596b7SLaurent Vivier /* Portions of this work are licensed under the terms of the GNU GPL, 18591596b7SLaurent Vivier * version 2 or later. See the COPYING file in the top-level directory. 19591596b7SLaurent Vivier */ 20591596b7SLaurent Vivier 21591596b7SLaurent Vivier #include "qemu/osdep.h" 22591596b7SLaurent Vivier #include "softfloat.h" 23591596b7SLaurent Vivier #include "fpu/softfloat-macros.h" 24591596b7SLaurent Vivier 25*0d379c17SLaurent Vivier static floatx80 propagateFloatx80NaNOneArg(floatx80 a, float_status *status) 26*0d379c17SLaurent Vivier { 27*0d379c17SLaurent Vivier if (floatx80_is_signaling_nan(a, status)) { 28*0d379c17SLaurent Vivier float_raise(float_flag_invalid, status); 29*0d379c17SLaurent Vivier } 30*0d379c17SLaurent Vivier 31*0d379c17SLaurent Vivier if (status->default_nan_mode) { 32*0d379c17SLaurent Vivier return floatx80_default_nan(status); 33*0d379c17SLaurent Vivier } 34*0d379c17SLaurent Vivier 35*0d379c17SLaurent Vivier return floatx80_maybe_silence_nan(a, status); 36*0d379c17SLaurent Vivier } 37*0d379c17SLaurent Vivier 38591596b7SLaurent Vivier /*---------------------------------------------------------------------------- 39591596b7SLaurent Vivier | Returns the modulo remainder of the extended double-precision floating-point 40591596b7SLaurent Vivier | value `a' with respect to the corresponding value `b'. 41591596b7SLaurent Vivier *----------------------------------------------------------------------------*/ 42591596b7SLaurent Vivier 43591596b7SLaurent Vivier floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status) 44591596b7SLaurent Vivier { 45591596b7SLaurent Vivier flag aSign, zSign; 46591596b7SLaurent Vivier int32_t aExp, bExp, expDiff; 47591596b7SLaurent Vivier uint64_t aSig0, aSig1, bSig; 48591596b7SLaurent Vivier uint64_t qTemp, term0, term1; 49591596b7SLaurent Vivier 50591596b7SLaurent Vivier aSig0 = extractFloatx80Frac(a); 51591596b7SLaurent Vivier aExp = extractFloatx80Exp(a); 52591596b7SLaurent Vivier aSign = extractFloatx80Sign(a); 53591596b7SLaurent Vivier bSig = extractFloatx80Frac(b); 54591596b7SLaurent Vivier bExp = extractFloatx80Exp(b); 55591596b7SLaurent Vivier 56591596b7SLaurent Vivier if (aExp == 0x7FFF) { 57591596b7SLaurent Vivier if ((uint64_t) (aSig0 << 1) 58591596b7SLaurent Vivier || ((bExp == 0x7FFF) && (uint64_t) (bSig << 1))) { 59591596b7SLaurent Vivier return propagateFloatx80NaN(a, b, status); 60591596b7SLaurent Vivier } 61591596b7SLaurent Vivier goto invalid; 62591596b7SLaurent Vivier } 63591596b7SLaurent Vivier if (bExp == 0x7FFF) { 64591596b7SLaurent Vivier if ((uint64_t) (bSig << 1)) { 65591596b7SLaurent Vivier return propagateFloatx80NaN(a, b, status); 66591596b7SLaurent Vivier } 67591596b7SLaurent Vivier return a; 68591596b7SLaurent Vivier } 69591596b7SLaurent Vivier if (bExp == 0) { 70591596b7SLaurent Vivier if (bSig == 0) { 71591596b7SLaurent Vivier invalid: 72591596b7SLaurent Vivier float_raise(float_flag_invalid, status); 73591596b7SLaurent Vivier return floatx80_default_nan(status); 74591596b7SLaurent Vivier } 75591596b7SLaurent Vivier normalizeFloatx80Subnormal(bSig, &bExp, &bSig); 76591596b7SLaurent Vivier } 77591596b7SLaurent Vivier if (aExp == 0) { 78591596b7SLaurent Vivier if ((uint64_t) (aSig0 << 1) == 0) { 79591596b7SLaurent Vivier return a; 80591596b7SLaurent Vivier } 81591596b7SLaurent Vivier normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0); 82591596b7SLaurent Vivier } 83591596b7SLaurent Vivier bSig |= LIT64(0x8000000000000000); 84591596b7SLaurent Vivier zSign = aSign; 85591596b7SLaurent Vivier expDiff = aExp - bExp; 86591596b7SLaurent Vivier aSig1 = 0; 87591596b7SLaurent Vivier if (expDiff < 0) { 88591596b7SLaurent Vivier return a; 89591596b7SLaurent Vivier } 90591596b7SLaurent Vivier qTemp = (bSig <= aSig0); 91591596b7SLaurent Vivier if (qTemp) { 92591596b7SLaurent Vivier aSig0 -= bSig; 93591596b7SLaurent Vivier } 94591596b7SLaurent Vivier expDiff -= 64; 95591596b7SLaurent Vivier while (0 < expDiff) { 96591596b7SLaurent Vivier qTemp = estimateDiv128To64(aSig0, aSig1, bSig); 97591596b7SLaurent Vivier qTemp = (2 < qTemp) ? qTemp - 2 : 0; 98591596b7SLaurent Vivier mul64To128(bSig, qTemp, &term0, &term1); 99591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 100591596b7SLaurent Vivier shortShift128Left(aSig0, aSig1, 62, &aSig0, &aSig1); 101591596b7SLaurent Vivier } 102591596b7SLaurent Vivier expDiff += 64; 103591596b7SLaurent Vivier if (0 < expDiff) { 104591596b7SLaurent Vivier qTemp = estimateDiv128To64(aSig0, aSig1, bSig); 105591596b7SLaurent Vivier qTemp = (2 < qTemp) ? qTemp - 2 : 0; 106591596b7SLaurent Vivier qTemp >>= 64 - expDiff; 107591596b7SLaurent Vivier mul64To128(bSig, qTemp << (64 - expDiff), &term0, &term1); 108591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 109591596b7SLaurent Vivier shortShift128Left(0, bSig, 64 - expDiff, &term0, &term1); 110591596b7SLaurent Vivier while (le128(term0, term1, aSig0, aSig1)) { 111591596b7SLaurent Vivier ++qTemp; 112591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 113591596b7SLaurent Vivier } 114591596b7SLaurent Vivier } 115591596b7SLaurent Vivier return 116591596b7SLaurent Vivier normalizeRoundAndPackFloatx80( 117591596b7SLaurent Vivier 80, zSign, bExp + expDiff, aSig0, aSig1, status); 118591596b7SLaurent Vivier } 119*0d379c17SLaurent Vivier 120*0d379c17SLaurent Vivier /*---------------------------------------------------------------------------- 121*0d379c17SLaurent Vivier | Returns the mantissa of the extended double-precision floating-point 122*0d379c17SLaurent Vivier | value `a'. 123*0d379c17SLaurent Vivier *----------------------------------------------------------------------------*/ 124*0d379c17SLaurent Vivier 125*0d379c17SLaurent Vivier floatx80 floatx80_getman(floatx80 a, float_status *status) 126*0d379c17SLaurent Vivier { 127*0d379c17SLaurent Vivier flag aSign; 128*0d379c17SLaurent Vivier int32_t aExp; 129*0d379c17SLaurent Vivier uint64_t aSig; 130*0d379c17SLaurent Vivier 131*0d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 132*0d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 133*0d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 134*0d379c17SLaurent Vivier 135*0d379c17SLaurent Vivier if (aExp == 0x7FFF) { 136*0d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 137*0d379c17SLaurent Vivier return propagateFloatx80NaNOneArg(a , status); 138*0d379c17SLaurent Vivier } 139*0d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 140*0d379c17SLaurent Vivier return floatx80_default_nan(status); 141*0d379c17SLaurent Vivier } 142*0d379c17SLaurent Vivier 143*0d379c17SLaurent Vivier if (aExp == 0) { 144*0d379c17SLaurent Vivier if (aSig == 0) { 145*0d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 146*0d379c17SLaurent Vivier } 147*0d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 148*0d379c17SLaurent Vivier } 149*0d379c17SLaurent Vivier 150*0d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign, 151*0d379c17SLaurent Vivier 0x3FFF, aSig, 0, status); 152*0d379c17SLaurent Vivier } 153*0d379c17SLaurent Vivier 154*0d379c17SLaurent Vivier /*---------------------------------------------------------------------------- 155*0d379c17SLaurent Vivier | Returns the exponent of the extended double-precision floating-point 156*0d379c17SLaurent Vivier | value `a' as an extended double-precision value. 157*0d379c17SLaurent Vivier *----------------------------------------------------------------------------*/ 158*0d379c17SLaurent Vivier 159*0d379c17SLaurent Vivier floatx80 floatx80_getexp(floatx80 a, float_status *status) 160*0d379c17SLaurent Vivier { 161*0d379c17SLaurent Vivier flag aSign; 162*0d379c17SLaurent Vivier int32_t aExp; 163*0d379c17SLaurent Vivier uint64_t aSig; 164*0d379c17SLaurent Vivier 165*0d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 166*0d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 167*0d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 168*0d379c17SLaurent Vivier 169*0d379c17SLaurent Vivier if (aExp == 0x7FFF) { 170*0d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 171*0d379c17SLaurent Vivier return propagateFloatx80NaNOneArg(a , status); 172*0d379c17SLaurent Vivier } 173*0d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 174*0d379c17SLaurent Vivier return floatx80_default_nan(status); 175*0d379c17SLaurent Vivier } 176*0d379c17SLaurent Vivier 177*0d379c17SLaurent Vivier if (aExp == 0) { 178*0d379c17SLaurent Vivier if (aSig == 0) { 179*0d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 180*0d379c17SLaurent Vivier } 181*0d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 182*0d379c17SLaurent Vivier } 183*0d379c17SLaurent Vivier 184*0d379c17SLaurent Vivier return int32_to_floatx80(aExp - 0x3FFF, status); 185*0d379c17SLaurent Vivier } 186*0d379c17SLaurent Vivier 187*0d379c17SLaurent Vivier /*---------------------------------------------------------------------------- 188*0d379c17SLaurent Vivier | Scales extended double-precision floating-point value in operand `a' by 189*0d379c17SLaurent Vivier | value `b'. The function truncates the value in the second operand 'b' to 190*0d379c17SLaurent Vivier | an integral value and adds that value to the exponent of the operand 'a'. 191*0d379c17SLaurent Vivier | The operation performed according to the IEC/IEEE Standard for Binary 192*0d379c17SLaurent Vivier | Floating-Point Arithmetic. 193*0d379c17SLaurent Vivier *----------------------------------------------------------------------------*/ 194*0d379c17SLaurent Vivier 195*0d379c17SLaurent Vivier floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status) 196*0d379c17SLaurent Vivier { 197*0d379c17SLaurent Vivier flag aSign, bSign; 198*0d379c17SLaurent Vivier int32_t aExp, bExp, shiftCount; 199*0d379c17SLaurent Vivier uint64_t aSig, bSig; 200*0d379c17SLaurent Vivier 201*0d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 202*0d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 203*0d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 204*0d379c17SLaurent Vivier bSig = extractFloatx80Frac(b); 205*0d379c17SLaurent Vivier bExp = extractFloatx80Exp(b); 206*0d379c17SLaurent Vivier bSign = extractFloatx80Sign(b); 207*0d379c17SLaurent Vivier 208*0d379c17SLaurent Vivier if (bExp == 0x7FFF) { 209*0d379c17SLaurent Vivier if ((uint64_t) (bSig << 1) || 210*0d379c17SLaurent Vivier ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) { 211*0d379c17SLaurent Vivier return propagateFloatx80NaN(a, b, status); 212*0d379c17SLaurent Vivier } 213*0d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 214*0d379c17SLaurent Vivier return floatx80_default_nan(status); 215*0d379c17SLaurent Vivier } 216*0d379c17SLaurent Vivier if (aExp == 0x7FFF) { 217*0d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 218*0d379c17SLaurent Vivier return propagateFloatx80NaN(a, b, status); 219*0d379c17SLaurent Vivier } 220*0d379c17SLaurent Vivier return packFloatx80(aSign, floatx80_infinity.high, 221*0d379c17SLaurent Vivier floatx80_infinity.low); 222*0d379c17SLaurent Vivier } 223*0d379c17SLaurent Vivier if (aExp == 0) { 224*0d379c17SLaurent Vivier if (aSig == 0) { 225*0d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 226*0d379c17SLaurent Vivier } 227*0d379c17SLaurent Vivier if (bExp < 0x3FFF) { 228*0d379c17SLaurent Vivier return a; 229*0d379c17SLaurent Vivier } 230*0d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 231*0d379c17SLaurent Vivier } 232*0d379c17SLaurent Vivier 233*0d379c17SLaurent Vivier if (bExp < 0x3FFF) { 234*0d379c17SLaurent Vivier return a; 235*0d379c17SLaurent Vivier } 236*0d379c17SLaurent Vivier 237*0d379c17SLaurent Vivier if (0x400F < bExp) { 238*0d379c17SLaurent Vivier aExp = bSign ? -0x6001 : 0xE000; 239*0d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 240*0d379c17SLaurent Vivier aSign, aExp, aSig, 0, status); 241*0d379c17SLaurent Vivier } 242*0d379c17SLaurent Vivier 243*0d379c17SLaurent Vivier shiftCount = 0x403E - bExp; 244*0d379c17SLaurent Vivier bSig >>= shiftCount; 245*0d379c17SLaurent Vivier aExp = bSign ? (aExp - bSig) : (aExp + bSig); 246*0d379c17SLaurent Vivier 247*0d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 248*0d379c17SLaurent Vivier aSign, aExp, aSig, 0, status); 249*0d379c17SLaurent Vivier } 250