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 17*808d77bcSLucien Murray-Pitts /* 18*808d77bcSLucien Murray-Pitts * Portions of this work are licensed under the terms of the GNU GPL, 19591596b7SLaurent Vivier * version 2 or later. See the COPYING file in the top-level directory. 20591596b7SLaurent Vivier */ 21591596b7SLaurent Vivier 22591596b7SLaurent Vivier #include "qemu/osdep.h" 23591596b7SLaurent Vivier #include "softfloat.h" 24591596b7SLaurent Vivier #include "fpu/softfloat-macros.h" 254b5c65b8SLaurent Vivier #include "softfloat_fpsp_tables.h" 26591596b7SLaurent Vivier 27c84813b8SLaurent Vivier #define pi_exp 0x4000 288c992abcSLaurent Vivier #define piby2_exp 0x3FFF 298c992abcSLaurent Vivier #define pi_sig LIT64(0xc90fdaa22168c235) 308c992abcSLaurent Vivier 310d379c17SLaurent Vivier static floatx80 propagateFloatx80NaNOneArg(floatx80 a, float_status *status) 320d379c17SLaurent Vivier { 330d379c17SLaurent Vivier if (floatx80_is_signaling_nan(a, status)) { 340d379c17SLaurent Vivier float_raise(float_flag_invalid, status); 351c0c951fSRichard Henderson a = floatx80_silence_nan(a, status); 360d379c17SLaurent Vivier } 370d379c17SLaurent Vivier 380d379c17SLaurent Vivier if (status->default_nan_mode) { 390d379c17SLaurent Vivier return floatx80_default_nan(status); 400d379c17SLaurent Vivier } 410d379c17SLaurent Vivier 421c0c951fSRichard Henderson return a; 430d379c17SLaurent Vivier } 440d379c17SLaurent Vivier 45*808d77bcSLucien Murray-Pitts /* 46*808d77bcSLucien Murray-Pitts * Returns the modulo remainder of the extended double-precision floating-point 47*808d77bcSLucien Murray-Pitts * value `a' with respect to the corresponding value `b'. 48*808d77bcSLucien Murray-Pitts */ 49591596b7SLaurent Vivier 50591596b7SLaurent Vivier floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status) 51591596b7SLaurent Vivier { 52591596b7SLaurent Vivier flag aSign, zSign; 53591596b7SLaurent Vivier int32_t aExp, bExp, expDiff; 54591596b7SLaurent Vivier uint64_t aSig0, aSig1, bSig; 55591596b7SLaurent Vivier uint64_t qTemp, term0, term1; 56591596b7SLaurent Vivier 57591596b7SLaurent Vivier aSig0 = extractFloatx80Frac(a); 58591596b7SLaurent Vivier aExp = extractFloatx80Exp(a); 59591596b7SLaurent Vivier aSign = extractFloatx80Sign(a); 60591596b7SLaurent Vivier bSig = extractFloatx80Frac(b); 61591596b7SLaurent Vivier bExp = extractFloatx80Exp(b); 62591596b7SLaurent Vivier 63591596b7SLaurent Vivier if (aExp == 0x7FFF) { 64591596b7SLaurent Vivier if ((uint64_t) (aSig0 << 1) 65591596b7SLaurent Vivier || ((bExp == 0x7FFF) && (uint64_t) (bSig << 1))) { 66591596b7SLaurent Vivier return propagateFloatx80NaN(a, b, status); 67591596b7SLaurent Vivier } 68591596b7SLaurent Vivier goto invalid; 69591596b7SLaurent Vivier } 70591596b7SLaurent Vivier if (bExp == 0x7FFF) { 71591596b7SLaurent Vivier if ((uint64_t) (bSig << 1)) { 72591596b7SLaurent Vivier return propagateFloatx80NaN(a, b, status); 73591596b7SLaurent Vivier } 74591596b7SLaurent Vivier return a; 75591596b7SLaurent Vivier } 76591596b7SLaurent Vivier if (bExp == 0) { 77591596b7SLaurent Vivier if (bSig == 0) { 78591596b7SLaurent Vivier invalid: 79591596b7SLaurent Vivier float_raise(float_flag_invalid, status); 80591596b7SLaurent Vivier return floatx80_default_nan(status); 81591596b7SLaurent Vivier } 82591596b7SLaurent Vivier normalizeFloatx80Subnormal(bSig, &bExp, &bSig); 83591596b7SLaurent Vivier } 84591596b7SLaurent Vivier if (aExp == 0) { 85591596b7SLaurent Vivier if ((uint64_t) (aSig0 << 1) == 0) { 86591596b7SLaurent Vivier return a; 87591596b7SLaurent Vivier } 88591596b7SLaurent Vivier normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0); 89591596b7SLaurent Vivier } 90591596b7SLaurent Vivier bSig |= LIT64(0x8000000000000000); 91591596b7SLaurent Vivier zSign = aSign; 92591596b7SLaurent Vivier expDiff = aExp - bExp; 93591596b7SLaurent Vivier aSig1 = 0; 94591596b7SLaurent Vivier if (expDiff < 0) { 95591596b7SLaurent Vivier return a; 96591596b7SLaurent Vivier } 97591596b7SLaurent Vivier qTemp = (bSig <= aSig0); 98591596b7SLaurent Vivier if (qTemp) { 99591596b7SLaurent Vivier aSig0 -= bSig; 100591596b7SLaurent Vivier } 101591596b7SLaurent Vivier expDiff -= 64; 102591596b7SLaurent Vivier while (0 < expDiff) { 103591596b7SLaurent Vivier qTemp = estimateDiv128To64(aSig0, aSig1, bSig); 104591596b7SLaurent Vivier qTemp = (2 < qTemp) ? qTemp - 2 : 0; 105591596b7SLaurent Vivier mul64To128(bSig, qTemp, &term0, &term1); 106591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 107591596b7SLaurent Vivier shortShift128Left(aSig0, aSig1, 62, &aSig0, &aSig1); 1085a73e7f3SLaurent Vivier expDiff -= 62; 109591596b7SLaurent Vivier } 110591596b7SLaurent Vivier expDiff += 64; 111591596b7SLaurent Vivier if (0 < expDiff) { 112591596b7SLaurent Vivier qTemp = estimateDiv128To64(aSig0, aSig1, bSig); 113591596b7SLaurent Vivier qTemp = (2 < qTemp) ? qTemp - 2 : 0; 114591596b7SLaurent Vivier qTemp >>= 64 - expDiff; 115591596b7SLaurent Vivier mul64To128(bSig, qTemp << (64 - expDiff), &term0, &term1); 116591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 117591596b7SLaurent Vivier shortShift128Left(0, bSig, 64 - expDiff, &term0, &term1); 118591596b7SLaurent Vivier while (le128(term0, term1, aSig0, aSig1)) { 119591596b7SLaurent Vivier ++qTemp; 120591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 121591596b7SLaurent Vivier } 122591596b7SLaurent Vivier } 123591596b7SLaurent Vivier return 124591596b7SLaurent Vivier normalizeRoundAndPackFloatx80( 125591596b7SLaurent Vivier 80, zSign, bExp + expDiff, aSig0, aSig1, status); 126591596b7SLaurent Vivier } 1270d379c17SLaurent Vivier 128*808d77bcSLucien Murray-Pitts /* 129*808d77bcSLucien Murray-Pitts * Returns the mantissa of the extended double-precision floating-point 130*808d77bcSLucien Murray-Pitts * value `a'. 131*808d77bcSLucien Murray-Pitts */ 1320d379c17SLaurent Vivier 1330d379c17SLaurent Vivier floatx80 floatx80_getman(floatx80 a, float_status *status) 1340d379c17SLaurent Vivier { 1350d379c17SLaurent Vivier flag aSign; 1360d379c17SLaurent Vivier int32_t aExp; 1370d379c17SLaurent Vivier uint64_t aSig; 1380d379c17SLaurent Vivier 1390d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 1400d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 1410d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 1420d379c17SLaurent Vivier 1430d379c17SLaurent Vivier if (aExp == 0x7FFF) { 1440d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 1450d379c17SLaurent Vivier return propagateFloatx80NaNOneArg(a , status); 1460d379c17SLaurent Vivier } 1470d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 1480d379c17SLaurent Vivier return floatx80_default_nan(status); 1490d379c17SLaurent Vivier } 1500d379c17SLaurent Vivier 1510d379c17SLaurent Vivier if (aExp == 0) { 1520d379c17SLaurent Vivier if (aSig == 0) { 1530d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 1540d379c17SLaurent Vivier } 1550d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 1560d379c17SLaurent Vivier } 1570d379c17SLaurent Vivier 1580d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign, 1590d379c17SLaurent Vivier 0x3FFF, aSig, 0, status); 1600d379c17SLaurent Vivier } 1610d379c17SLaurent Vivier 162*808d77bcSLucien Murray-Pitts /* 163*808d77bcSLucien Murray-Pitts * Returns the exponent of the extended double-precision floating-point 164*808d77bcSLucien Murray-Pitts * value `a' as an extended double-precision value. 165*808d77bcSLucien Murray-Pitts */ 1660d379c17SLaurent Vivier 1670d379c17SLaurent Vivier floatx80 floatx80_getexp(floatx80 a, float_status *status) 1680d379c17SLaurent Vivier { 1690d379c17SLaurent Vivier flag aSign; 1700d379c17SLaurent Vivier int32_t aExp; 1710d379c17SLaurent Vivier uint64_t aSig; 1720d379c17SLaurent Vivier 1730d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 1740d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 1750d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 1760d379c17SLaurent Vivier 1770d379c17SLaurent Vivier if (aExp == 0x7FFF) { 1780d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 1790d379c17SLaurent Vivier return propagateFloatx80NaNOneArg(a , status); 1800d379c17SLaurent Vivier } 1810d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 1820d379c17SLaurent Vivier return floatx80_default_nan(status); 1830d379c17SLaurent Vivier } 1840d379c17SLaurent Vivier 1850d379c17SLaurent Vivier if (aExp == 0) { 1860d379c17SLaurent Vivier if (aSig == 0) { 1870d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 1880d379c17SLaurent Vivier } 1890d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 1900d379c17SLaurent Vivier } 1910d379c17SLaurent Vivier 1920d379c17SLaurent Vivier return int32_to_floatx80(aExp - 0x3FFF, status); 1930d379c17SLaurent Vivier } 1940d379c17SLaurent Vivier 195*808d77bcSLucien Murray-Pitts /* 196*808d77bcSLucien Murray-Pitts * Scales extended double-precision floating-point value in operand `a' by 197*808d77bcSLucien Murray-Pitts * value `b'. The function truncates the value in the second operand 'b' to 198*808d77bcSLucien Murray-Pitts * an integral value and adds that value to the exponent of the operand 'a'. 199*808d77bcSLucien Murray-Pitts * The operation performed according to the IEC/IEEE Standard for Binary 200*808d77bcSLucien Murray-Pitts * Floating-Point Arithmetic. 201*808d77bcSLucien Murray-Pitts */ 2020d379c17SLaurent Vivier 2030d379c17SLaurent Vivier floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status) 2040d379c17SLaurent Vivier { 2050d379c17SLaurent Vivier flag aSign, bSign; 2060d379c17SLaurent Vivier int32_t aExp, bExp, shiftCount; 2070d379c17SLaurent Vivier uint64_t aSig, bSig; 2080d379c17SLaurent Vivier 2090d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 2100d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 2110d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 2120d379c17SLaurent Vivier bSig = extractFloatx80Frac(b); 2130d379c17SLaurent Vivier bExp = extractFloatx80Exp(b); 2140d379c17SLaurent Vivier bSign = extractFloatx80Sign(b); 2150d379c17SLaurent Vivier 2160d379c17SLaurent Vivier if (bExp == 0x7FFF) { 2170d379c17SLaurent Vivier if ((uint64_t) (bSig << 1) || 2180d379c17SLaurent Vivier ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) { 2190d379c17SLaurent Vivier return propagateFloatx80NaN(a, b, status); 2200d379c17SLaurent Vivier } 2210d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 2220d379c17SLaurent Vivier return floatx80_default_nan(status); 2230d379c17SLaurent Vivier } 2240d379c17SLaurent Vivier if (aExp == 0x7FFF) { 2250d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 2260d379c17SLaurent Vivier return propagateFloatx80NaN(a, b, status); 2270d379c17SLaurent Vivier } 2280d379c17SLaurent Vivier return packFloatx80(aSign, floatx80_infinity.high, 2290d379c17SLaurent Vivier floatx80_infinity.low); 2300d379c17SLaurent Vivier } 2310d379c17SLaurent Vivier if (aExp == 0) { 2320d379c17SLaurent Vivier if (aSig == 0) { 2330d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 2340d379c17SLaurent Vivier } 2350d379c17SLaurent Vivier if (bExp < 0x3FFF) { 2360d379c17SLaurent Vivier return a; 2370d379c17SLaurent Vivier } 2380d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 2390d379c17SLaurent Vivier } 2400d379c17SLaurent Vivier 2410d379c17SLaurent Vivier if (bExp < 0x3FFF) { 2420d379c17SLaurent Vivier return a; 2430d379c17SLaurent Vivier } 2440d379c17SLaurent Vivier 2450d379c17SLaurent Vivier if (0x400F < bExp) { 2460d379c17SLaurent Vivier aExp = bSign ? -0x6001 : 0xE000; 2470d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 2480d379c17SLaurent Vivier aSign, aExp, aSig, 0, status); 2490d379c17SLaurent Vivier } 2500d379c17SLaurent Vivier 2510d379c17SLaurent Vivier shiftCount = 0x403E - bExp; 2520d379c17SLaurent Vivier bSig >>= shiftCount; 2530d379c17SLaurent Vivier aExp = bSign ? (aExp - bSig) : (aExp + bSig); 2540d379c17SLaurent Vivier 2550d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 2560d379c17SLaurent Vivier aSign, aExp, aSig, 0, status); 2570d379c17SLaurent Vivier } 2589a069775SLaurent Vivier 2599a069775SLaurent Vivier floatx80 floatx80_move(floatx80 a, float_status *status) 2609a069775SLaurent Vivier { 2619a069775SLaurent Vivier flag aSign; 2629a069775SLaurent Vivier int32_t aExp; 2639a069775SLaurent Vivier uint64_t aSig; 2649a069775SLaurent Vivier 2659a069775SLaurent Vivier aSig = extractFloatx80Frac(a); 2669a069775SLaurent Vivier aExp = extractFloatx80Exp(a); 2679a069775SLaurent Vivier aSign = extractFloatx80Sign(a); 2689a069775SLaurent Vivier 2699a069775SLaurent Vivier if (aExp == 0x7FFF) { 2709a069775SLaurent Vivier if ((uint64_t)(aSig << 1)) { 2719a069775SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 2729a069775SLaurent Vivier } 2739a069775SLaurent Vivier return a; 2749a069775SLaurent Vivier } 2759a069775SLaurent Vivier if (aExp == 0) { 2769a069775SLaurent Vivier if (aSig == 0) { 2779a069775SLaurent Vivier return a; 2789a069775SLaurent Vivier } 2799a069775SLaurent Vivier normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision, 2809a069775SLaurent Vivier aSign, aExp, aSig, 0, status); 2819a069775SLaurent Vivier } 2829a069775SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign, 2839a069775SLaurent Vivier aExp, aSig, 0, status); 2849a069775SLaurent Vivier } 2854b5c65b8SLaurent Vivier 286*808d77bcSLucien Murray-Pitts /* 287*808d77bcSLucien Murray-Pitts * Algorithms for transcendental functions supported by MC68881 and MC68882 288*808d77bcSLucien Murray-Pitts * mathematical coprocessors. The functions are derived from FPSP library. 289*808d77bcSLucien Murray-Pitts */ 2904b5c65b8SLaurent Vivier 2914b5c65b8SLaurent Vivier #define one_exp 0x3FFF 2924b5c65b8SLaurent Vivier #define one_sig LIT64(0x8000000000000000) 2934b5c65b8SLaurent Vivier 294*808d77bcSLucien Murray-Pitts /* 295*808d77bcSLucien Murray-Pitts * Function for compactifying extended double-precision floating point values. 296*808d77bcSLucien Murray-Pitts */ 2974b5c65b8SLaurent Vivier 2984b5c65b8SLaurent Vivier static int32_t floatx80_make_compact(int32_t aExp, uint64_t aSig) 2994b5c65b8SLaurent Vivier { 3004b5c65b8SLaurent Vivier return (aExp << 16) | (aSig >> 48); 3014b5c65b8SLaurent Vivier } 3024b5c65b8SLaurent Vivier 303*808d77bcSLucien Murray-Pitts /* 304*808d77bcSLucien Murray-Pitts * Log base e of x plus 1 305*808d77bcSLucien Murray-Pitts */ 3064b5c65b8SLaurent Vivier 3074b5c65b8SLaurent Vivier floatx80 floatx80_lognp1(floatx80 a, float_status *status) 3084b5c65b8SLaurent Vivier { 3094b5c65b8SLaurent Vivier flag aSign; 3104b5c65b8SLaurent Vivier int32_t aExp; 3114b5c65b8SLaurent Vivier uint64_t aSig, fSig; 3124b5c65b8SLaurent Vivier 3134b5c65b8SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 3144b5c65b8SLaurent Vivier 3154b5c65b8SLaurent Vivier int32_t compact, j, k; 3164b5c65b8SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu; 3174b5c65b8SLaurent Vivier 3184b5c65b8SLaurent Vivier aSig = extractFloatx80Frac(a); 3194b5c65b8SLaurent Vivier aExp = extractFloatx80Exp(a); 3204b5c65b8SLaurent Vivier aSign = extractFloatx80Sign(a); 3214b5c65b8SLaurent Vivier 3224b5c65b8SLaurent Vivier if (aExp == 0x7FFF) { 3234b5c65b8SLaurent Vivier if ((uint64_t) (aSig << 1)) { 3244b5c65b8SLaurent Vivier propagateFloatx80NaNOneArg(a, status); 3254b5c65b8SLaurent Vivier } 3264b5c65b8SLaurent Vivier if (aSign) { 3274b5c65b8SLaurent Vivier float_raise(float_flag_invalid, status); 3284b5c65b8SLaurent Vivier return floatx80_default_nan(status); 3294b5c65b8SLaurent Vivier } 3304b5c65b8SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, floatx80_infinity.low); 3314b5c65b8SLaurent Vivier } 3324b5c65b8SLaurent Vivier 3334b5c65b8SLaurent Vivier if (aExp == 0 && aSig == 0) { 3344b5c65b8SLaurent Vivier return packFloatx80(aSign, 0, 0); 3354b5c65b8SLaurent Vivier } 3364b5c65b8SLaurent Vivier 3374b5c65b8SLaurent Vivier if (aSign && aExp >= one_exp) { 3384b5c65b8SLaurent Vivier if (aExp == one_exp && aSig == one_sig) { 3394b5c65b8SLaurent Vivier float_raise(float_flag_divbyzero, status); 340981348afSLaurent Vivier return packFloatx80(aSign, floatx80_infinity.high, 341981348afSLaurent Vivier floatx80_infinity.low); 3424b5c65b8SLaurent Vivier } 3434b5c65b8SLaurent Vivier float_raise(float_flag_invalid, status); 3444b5c65b8SLaurent Vivier return floatx80_default_nan(status); 3454b5c65b8SLaurent Vivier } 3464b5c65b8SLaurent Vivier 3474b5c65b8SLaurent Vivier if (aExp < 0x3f99 || (aExp == 0x3f99 && aSig == one_sig)) { 3484b5c65b8SLaurent Vivier /* <= min threshold */ 3494b5c65b8SLaurent Vivier float_raise(float_flag_inexact, status); 3504b5c65b8SLaurent Vivier return floatx80_move(a, status); 3514b5c65b8SLaurent Vivier } 3524b5c65b8SLaurent Vivier 3534b5c65b8SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 3544b5c65b8SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 3554b5c65b8SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 3564b5c65b8SLaurent Vivier status->floatx80_rounding_precision = 80; 3574b5c65b8SLaurent Vivier 3584b5c65b8SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 3594b5c65b8SLaurent Vivier 3604b5c65b8SLaurent Vivier fp0 = a; /* Z */ 3614b5c65b8SLaurent Vivier fp1 = a; 3624b5c65b8SLaurent Vivier 3634b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000), 3644b5c65b8SLaurent Vivier status), status); /* X = (1+Z) */ 3654b5c65b8SLaurent Vivier 3664b5c65b8SLaurent Vivier aExp = extractFloatx80Exp(fp0); 3674b5c65b8SLaurent Vivier aSig = extractFloatx80Frac(fp0); 3684b5c65b8SLaurent Vivier 3694b5c65b8SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 3704b5c65b8SLaurent Vivier 3714b5c65b8SLaurent Vivier if (compact < 0x3FFE8000 || compact > 0x3FFFC000) { 3724b5c65b8SLaurent Vivier /* |X| < 1/2 or |X| > 3/2 */ 3734b5c65b8SLaurent Vivier k = aExp - 0x3FFF; 3744b5c65b8SLaurent Vivier fp1 = int32_to_floatx80(k, status); 3754b5c65b8SLaurent Vivier 3764b5c65b8SLaurent Vivier fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); 3774b5c65b8SLaurent Vivier j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */ 3784b5c65b8SLaurent Vivier 3794b5c65b8SLaurent Vivier f = packFloatx80(0, 0x3FFF, fSig); /* F */ 3804b5c65b8SLaurent Vivier fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */ 3814b5c65b8SLaurent Vivier 3824b5c65b8SLaurent Vivier fp0 = floatx80_sub(fp0, f, status); /* Y-F */ 3834b5c65b8SLaurent Vivier 3844b5c65b8SLaurent Vivier lp1cont1: 3854b5c65b8SLaurent Vivier /* LP1CONT1 */ 3864b5c65b8SLaurent Vivier fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */ 3874b5c65b8SLaurent Vivier logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); 3884b5c65b8SLaurent Vivier klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */ 3894b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */ 3904b5c65b8SLaurent Vivier 3914b5c65b8SLaurent Vivier fp3 = fp2; 3924b5c65b8SLaurent Vivier fp1 = fp2; 3934b5c65b8SLaurent Vivier 3944b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 3954b5c65b8SLaurent Vivier make_float64(0x3FC2499AB5E4040B), status), 3964b5c65b8SLaurent Vivier status); /* V*A6 */ 3974b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, float64_to_floatx80( 3984b5c65b8SLaurent Vivier make_float64(0xBFC555B5848CB7DB), status), 3994b5c65b8SLaurent Vivier status); /* V*A5 */ 4004b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 4014b5c65b8SLaurent Vivier make_float64(0x3FC99999987D8730), status), 4024b5c65b8SLaurent Vivier status); /* A4+V*A6 */ 4034b5c65b8SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 4044b5c65b8SLaurent Vivier make_float64(0xBFCFFFFFFF6F7E97), status), 4054b5c65b8SLaurent Vivier status); /* A3+V*A5 */ 4064b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */ 4074b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */ 4084b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 4094b5c65b8SLaurent Vivier make_float64(0x3FD55555555555A4), status), 4104b5c65b8SLaurent Vivier status); /* A2+V*(A4+V*A6) */ 4114b5c65b8SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 4124b5c65b8SLaurent Vivier make_float64(0xBFE0000000000008), status), 4134b5c65b8SLaurent Vivier status); /* A1+V*(A3+V*A5) */ 4144b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */ 4154b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */ 4164b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */ 4174b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */ 4184b5c65b8SLaurent Vivier 4194b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, log_tbl[j + 1], 4204b5c65b8SLaurent Vivier status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */ 4214b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */ 4224b5c65b8SLaurent Vivier 4234b5c65b8SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 4244b5c65b8SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 4254b5c65b8SLaurent Vivier 4264b5c65b8SLaurent Vivier a = floatx80_add(fp0, klog2, status); 4274b5c65b8SLaurent Vivier 4284b5c65b8SLaurent Vivier float_raise(float_flag_inexact, status); 4294b5c65b8SLaurent Vivier 4304b5c65b8SLaurent Vivier return a; 4314b5c65b8SLaurent Vivier } else if (compact < 0x3FFEF07D || compact > 0x3FFF8841) { 4324b5c65b8SLaurent Vivier /* |X| < 1/16 or |X| > -1/16 */ 4334b5c65b8SLaurent Vivier /* LP1CARE */ 4344b5c65b8SLaurent Vivier fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); 4354b5c65b8SLaurent Vivier f = packFloatx80(0, 0x3FFF, fSig); /* F */ 4364b5c65b8SLaurent Vivier j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */ 4374b5c65b8SLaurent Vivier 4384b5c65b8SLaurent Vivier if (compact >= 0x3FFF8000) { /* 1+Z >= 1 */ 4394b5c65b8SLaurent Vivier /* KISZERO */ 4404b5c65b8SLaurent Vivier fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x3F800000), 4414b5c65b8SLaurent Vivier status), f, status); /* 1-F */ 4424b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (1-F)+Z */ 4434b5c65b8SLaurent Vivier fp1 = packFloatx80(0, 0, 0); /* K = 0 */ 4444b5c65b8SLaurent Vivier } else { 4454b5c65b8SLaurent Vivier /* KISNEG */ 4464b5c65b8SLaurent Vivier fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x40000000), 4474b5c65b8SLaurent Vivier status), f, status); /* 2-F */ 4484b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, fp1, status); /* 2Z */ 4494b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (2-F)+2Z */ 4504b5c65b8SLaurent Vivier fp1 = packFloatx80(1, one_exp, one_sig); /* K = -1 */ 4514b5c65b8SLaurent Vivier } 4524b5c65b8SLaurent Vivier goto lp1cont1; 4534b5c65b8SLaurent Vivier } else { 4544b5c65b8SLaurent Vivier /* LP1ONE16 */ 4554b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2Z */ 4564b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000), 4574b5c65b8SLaurent Vivier status), status); /* FP0 IS 1+X */ 4584b5c65b8SLaurent Vivier 4594b5c65b8SLaurent Vivier /* LP1CONT2 */ 4604b5c65b8SLaurent Vivier fp1 = floatx80_div(fp1, fp0, status); /* U */ 4614b5c65b8SLaurent Vivier saveu = fp1; 4624b5c65b8SLaurent Vivier fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */ 4634b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */ 4644b5c65b8SLaurent Vivier 4654b5c65b8SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6), 4664b5c65b8SLaurent Vivier status); /* B5 */ 4674b5c65b8SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0), 4684b5c65b8SLaurent Vivier status); /* B4 */ 4694b5c65b8SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */ 4704b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */ 4714b5c65b8SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 4724b5c65b8SLaurent Vivier make_float64(0x3F624924928BCCFF), status), 4734b5c65b8SLaurent Vivier status); /* B3+W*B5 */ 4744b5c65b8SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 4754b5c65b8SLaurent Vivier make_float64(0x3F899999999995EC), status), 4764b5c65b8SLaurent Vivier status); /* B2+W*B4 */ 4774b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */ 4784b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */ 4794b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 4804b5c65b8SLaurent Vivier make_float64(0x3FB5555555555555), status), 4814b5c65b8SLaurent Vivier status); /* B1+W*(B3+W*B5) */ 4824b5c65b8SLaurent Vivier 4834b5c65b8SLaurent Vivier fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */ 4844b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, fp2, 4854b5c65b8SLaurent Vivier status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */ 4864b5c65b8SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, 4874b5c65b8SLaurent Vivier status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */ 4884b5c65b8SLaurent Vivier 4894b5c65b8SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 4904b5c65b8SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 4914b5c65b8SLaurent Vivier 4924b5c65b8SLaurent Vivier a = floatx80_add(fp0, saveu, status); 4934b5c65b8SLaurent Vivier 4944b5c65b8SLaurent Vivier /*if (!floatx80_is_zero(a)) { */ 4954b5c65b8SLaurent Vivier float_raise(float_flag_inexact, status); 4964b5c65b8SLaurent Vivier /*} */ 4974b5c65b8SLaurent Vivier 4984b5c65b8SLaurent Vivier return a; 4994b5c65b8SLaurent Vivier } 5004b5c65b8SLaurent Vivier } 50150067bd1SLaurent Vivier 502*808d77bcSLucien Murray-Pitts /* 503*808d77bcSLucien Murray-Pitts * Log base e 504*808d77bcSLucien Murray-Pitts */ 50550067bd1SLaurent Vivier 50650067bd1SLaurent Vivier floatx80 floatx80_logn(floatx80 a, float_status *status) 50750067bd1SLaurent Vivier { 50850067bd1SLaurent Vivier flag aSign; 50950067bd1SLaurent Vivier int32_t aExp; 51050067bd1SLaurent Vivier uint64_t aSig, fSig; 51150067bd1SLaurent Vivier 51250067bd1SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 51350067bd1SLaurent Vivier 51450067bd1SLaurent Vivier int32_t compact, j, k, adjk; 51550067bd1SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu; 51650067bd1SLaurent Vivier 51750067bd1SLaurent Vivier aSig = extractFloatx80Frac(a); 51850067bd1SLaurent Vivier aExp = extractFloatx80Exp(a); 51950067bd1SLaurent Vivier aSign = extractFloatx80Sign(a); 52050067bd1SLaurent Vivier 52150067bd1SLaurent Vivier if (aExp == 0x7FFF) { 52250067bd1SLaurent Vivier if ((uint64_t) (aSig << 1)) { 52350067bd1SLaurent Vivier propagateFloatx80NaNOneArg(a, status); 52450067bd1SLaurent Vivier } 52550067bd1SLaurent Vivier if (aSign == 0) { 52650067bd1SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 52750067bd1SLaurent Vivier floatx80_infinity.low); 52850067bd1SLaurent Vivier } 52950067bd1SLaurent Vivier } 53050067bd1SLaurent Vivier 53150067bd1SLaurent Vivier adjk = 0; 53250067bd1SLaurent Vivier 53350067bd1SLaurent Vivier if (aExp == 0) { 53450067bd1SLaurent Vivier if (aSig == 0) { /* zero */ 53550067bd1SLaurent Vivier float_raise(float_flag_divbyzero, status); 53650067bd1SLaurent Vivier return packFloatx80(1, floatx80_infinity.high, 53750067bd1SLaurent Vivier floatx80_infinity.low); 53850067bd1SLaurent Vivier } 53950067bd1SLaurent Vivier if ((aSig & one_sig) == 0) { /* denormal */ 54050067bd1SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 54150067bd1SLaurent Vivier adjk = -100; 54250067bd1SLaurent Vivier aExp += 100; 54350067bd1SLaurent Vivier a = packFloatx80(aSign, aExp, aSig); 54450067bd1SLaurent Vivier } 54550067bd1SLaurent Vivier } 54650067bd1SLaurent Vivier 54750067bd1SLaurent Vivier if (aSign) { 54850067bd1SLaurent Vivier float_raise(float_flag_invalid, status); 54950067bd1SLaurent Vivier return floatx80_default_nan(status); 55050067bd1SLaurent Vivier } 55150067bd1SLaurent Vivier 55250067bd1SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 55350067bd1SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 55450067bd1SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 55550067bd1SLaurent Vivier status->floatx80_rounding_precision = 80; 55650067bd1SLaurent Vivier 55750067bd1SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 55850067bd1SLaurent Vivier 55950067bd1SLaurent Vivier if (compact < 0x3FFEF07D || compact > 0x3FFF8841) { 56050067bd1SLaurent Vivier /* |X| < 15/16 or |X| > 17/16 */ 56150067bd1SLaurent Vivier k = aExp - 0x3FFF; 56250067bd1SLaurent Vivier k += adjk; 56350067bd1SLaurent Vivier fp1 = int32_to_floatx80(k, status); 56450067bd1SLaurent Vivier 56550067bd1SLaurent Vivier fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); 56650067bd1SLaurent Vivier j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */ 56750067bd1SLaurent Vivier 56850067bd1SLaurent Vivier f = packFloatx80(0, 0x3FFF, fSig); /* F */ 56950067bd1SLaurent Vivier fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */ 57050067bd1SLaurent Vivier 57150067bd1SLaurent Vivier fp0 = floatx80_sub(fp0, f, status); /* Y-F */ 57250067bd1SLaurent Vivier 57350067bd1SLaurent Vivier /* LP1CONT1 */ 57450067bd1SLaurent Vivier fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */ 57550067bd1SLaurent Vivier logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); 57650067bd1SLaurent Vivier klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */ 57750067bd1SLaurent Vivier fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */ 57850067bd1SLaurent Vivier 57950067bd1SLaurent Vivier fp3 = fp2; 58050067bd1SLaurent Vivier fp1 = fp2; 58150067bd1SLaurent Vivier 58250067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 58350067bd1SLaurent Vivier make_float64(0x3FC2499AB5E4040B), status), 58450067bd1SLaurent Vivier status); /* V*A6 */ 58550067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, float64_to_floatx80( 58650067bd1SLaurent Vivier make_float64(0xBFC555B5848CB7DB), status), 58750067bd1SLaurent Vivier status); /* V*A5 */ 58850067bd1SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 58950067bd1SLaurent Vivier make_float64(0x3FC99999987D8730), status), 59050067bd1SLaurent Vivier status); /* A4+V*A6 */ 59150067bd1SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 59250067bd1SLaurent Vivier make_float64(0xBFCFFFFFFF6F7E97), status), 59350067bd1SLaurent Vivier status); /* A3+V*A5 */ 59450067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */ 59550067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */ 59650067bd1SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 59750067bd1SLaurent Vivier make_float64(0x3FD55555555555A4), status), 59850067bd1SLaurent Vivier status); /* A2+V*(A4+V*A6) */ 59950067bd1SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 60050067bd1SLaurent Vivier make_float64(0xBFE0000000000008), status), 60150067bd1SLaurent Vivier status); /* A1+V*(A3+V*A5) */ 60250067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */ 60350067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */ 60450067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */ 60550067bd1SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */ 60650067bd1SLaurent Vivier 60750067bd1SLaurent Vivier fp1 = floatx80_add(fp1, log_tbl[j + 1], 60850067bd1SLaurent Vivier status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */ 60950067bd1SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */ 61050067bd1SLaurent Vivier 61150067bd1SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 61250067bd1SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 61350067bd1SLaurent Vivier 61450067bd1SLaurent Vivier a = floatx80_add(fp0, klog2, status); 61550067bd1SLaurent Vivier 61650067bd1SLaurent Vivier float_raise(float_flag_inexact, status); 61750067bd1SLaurent Vivier 61850067bd1SLaurent Vivier return a; 61950067bd1SLaurent Vivier } else { /* |X-1| >= 1/16 */ 62050067bd1SLaurent Vivier fp0 = a; 62150067bd1SLaurent Vivier fp1 = a; 62250067bd1SLaurent Vivier fp1 = floatx80_sub(fp1, float32_to_floatx80(make_float32(0x3F800000), 62350067bd1SLaurent Vivier status), status); /* FP1 IS X-1 */ 62450067bd1SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000), 62550067bd1SLaurent Vivier status), status); /* FP0 IS X+1 */ 62650067bd1SLaurent Vivier fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2(X-1) */ 62750067bd1SLaurent Vivier 62850067bd1SLaurent Vivier /* LP1CONT2 */ 62950067bd1SLaurent Vivier fp1 = floatx80_div(fp1, fp0, status); /* U */ 63050067bd1SLaurent Vivier saveu = fp1; 63150067bd1SLaurent Vivier fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */ 63250067bd1SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */ 63350067bd1SLaurent Vivier 63450067bd1SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6), 63550067bd1SLaurent Vivier status); /* B5 */ 63650067bd1SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0), 63750067bd1SLaurent Vivier status); /* B4 */ 63850067bd1SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */ 63950067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */ 64050067bd1SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 64150067bd1SLaurent Vivier make_float64(0x3F624924928BCCFF), status), 64250067bd1SLaurent Vivier status); /* B3+W*B5 */ 64350067bd1SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 64450067bd1SLaurent Vivier make_float64(0x3F899999999995EC), status), 64550067bd1SLaurent Vivier status); /* B2+W*B4 */ 64650067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */ 64750067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */ 64850067bd1SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 64950067bd1SLaurent Vivier make_float64(0x3FB5555555555555), status), 65050067bd1SLaurent Vivier status); /* B1+W*(B3+W*B5) */ 65150067bd1SLaurent Vivier 65250067bd1SLaurent Vivier fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */ 65350067bd1SLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */ 65450067bd1SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, 65550067bd1SLaurent Vivier status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */ 65650067bd1SLaurent Vivier 65750067bd1SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 65850067bd1SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 65950067bd1SLaurent Vivier 66050067bd1SLaurent Vivier a = floatx80_add(fp0, saveu, status); 66150067bd1SLaurent Vivier 66250067bd1SLaurent Vivier /*if (!floatx80_is_zero(a)) { */ 66350067bd1SLaurent Vivier float_raise(float_flag_inexact, status); 66450067bd1SLaurent Vivier /*} */ 66550067bd1SLaurent Vivier 66650067bd1SLaurent Vivier return a; 66750067bd1SLaurent Vivier } 66850067bd1SLaurent Vivier } 669248efb66SLaurent Vivier 670*808d77bcSLucien Murray-Pitts /* 671*808d77bcSLucien Murray-Pitts * Log base 10 672*808d77bcSLucien Murray-Pitts */ 673248efb66SLaurent Vivier 674248efb66SLaurent Vivier floatx80 floatx80_log10(floatx80 a, float_status *status) 675248efb66SLaurent Vivier { 676248efb66SLaurent Vivier flag aSign; 677248efb66SLaurent Vivier int32_t aExp; 678248efb66SLaurent Vivier uint64_t aSig; 679248efb66SLaurent Vivier 680248efb66SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 681248efb66SLaurent Vivier 682248efb66SLaurent Vivier floatx80 fp0, fp1; 683248efb66SLaurent Vivier 684248efb66SLaurent Vivier aSig = extractFloatx80Frac(a); 685248efb66SLaurent Vivier aExp = extractFloatx80Exp(a); 686248efb66SLaurent Vivier aSign = extractFloatx80Sign(a); 687248efb66SLaurent Vivier 688248efb66SLaurent Vivier if (aExp == 0x7FFF) { 689248efb66SLaurent Vivier if ((uint64_t) (aSig << 1)) { 690248efb66SLaurent Vivier propagateFloatx80NaNOneArg(a, status); 691248efb66SLaurent Vivier } 692248efb66SLaurent Vivier if (aSign == 0) { 693248efb66SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 694248efb66SLaurent Vivier floatx80_infinity.low); 695248efb66SLaurent Vivier } 696248efb66SLaurent Vivier } 697248efb66SLaurent Vivier 698248efb66SLaurent Vivier if (aExp == 0 && aSig == 0) { 699248efb66SLaurent Vivier float_raise(float_flag_divbyzero, status); 700248efb66SLaurent Vivier return packFloatx80(1, floatx80_infinity.high, 701248efb66SLaurent Vivier floatx80_infinity.low); 702248efb66SLaurent Vivier } 703248efb66SLaurent Vivier 704248efb66SLaurent Vivier if (aSign) { 705248efb66SLaurent Vivier float_raise(float_flag_invalid, status); 706248efb66SLaurent Vivier return floatx80_default_nan(status); 707248efb66SLaurent Vivier } 708248efb66SLaurent Vivier 709248efb66SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 710248efb66SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 711248efb66SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 712248efb66SLaurent Vivier status->floatx80_rounding_precision = 80; 713248efb66SLaurent Vivier 714248efb66SLaurent Vivier fp0 = floatx80_logn(a, status); 715248efb66SLaurent Vivier fp1 = packFloatx80(0, 0x3FFD, LIT64(0xDE5BD8A937287195)); /* INV_L10 */ 716248efb66SLaurent Vivier 717248efb66SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 718248efb66SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 719248efb66SLaurent Vivier 720248efb66SLaurent Vivier a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L10 */ 721248efb66SLaurent Vivier 722248efb66SLaurent Vivier float_raise(float_flag_inexact, status); 723248efb66SLaurent Vivier 724248efb66SLaurent Vivier return a; 725248efb66SLaurent Vivier } 72667b453edSLaurent Vivier 727*808d77bcSLucien Murray-Pitts /* 728*808d77bcSLucien Murray-Pitts * Log base 2 729*808d77bcSLucien Murray-Pitts */ 73067b453edSLaurent Vivier 73167b453edSLaurent Vivier floatx80 floatx80_log2(floatx80 a, float_status *status) 73267b453edSLaurent Vivier { 73367b453edSLaurent Vivier flag aSign; 73467b453edSLaurent Vivier int32_t aExp; 73567b453edSLaurent Vivier uint64_t aSig; 73667b453edSLaurent Vivier 73767b453edSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 73867b453edSLaurent Vivier 73967b453edSLaurent Vivier floatx80 fp0, fp1; 74067b453edSLaurent Vivier 74167b453edSLaurent Vivier aSig = extractFloatx80Frac(a); 74267b453edSLaurent Vivier aExp = extractFloatx80Exp(a); 74367b453edSLaurent Vivier aSign = extractFloatx80Sign(a); 74467b453edSLaurent Vivier 74567b453edSLaurent Vivier if (aExp == 0x7FFF) { 74667b453edSLaurent Vivier if ((uint64_t) (aSig << 1)) { 74767b453edSLaurent Vivier propagateFloatx80NaNOneArg(a, status); 74867b453edSLaurent Vivier } 74967b453edSLaurent Vivier if (aSign == 0) { 75067b453edSLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 75167b453edSLaurent Vivier floatx80_infinity.low); 75267b453edSLaurent Vivier } 75367b453edSLaurent Vivier } 75467b453edSLaurent Vivier 75567b453edSLaurent Vivier if (aExp == 0) { 75667b453edSLaurent Vivier if (aSig == 0) { 75767b453edSLaurent Vivier float_raise(float_flag_divbyzero, status); 75867b453edSLaurent Vivier return packFloatx80(1, floatx80_infinity.high, 75967b453edSLaurent Vivier floatx80_infinity.low); 76067b453edSLaurent Vivier } 76167b453edSLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 76267b453edSLaurent Vivier } 76367b453edSLaurent Vivier 76467b453edSLaurent Vivier if (aSign) { 76567b453edSLaurent Vivier float_raise(float_flag_invalid, status); 76667b453edSLaurent Vivier return floatx80_default_nan(status); 76767b453edSLaurent Vivier } 76867b453edSLaurent Vivier 76967b453edSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 77067b453edSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 77167b453edSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 77267b453edSLaurent Vivier status->floatx80_rounding_precision = 80; 77367b453edSLaurent Vivier 77467b453edSLaurent Vivier if (aSig == one_sig) { /* X is 2^k */ 77567b453edSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 77667b453edSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 77767b453edSLaurent Vivier 77867b453edSLaurent Vivier a = int32_to_floatx80(aExp - 0x3FFF, status); 77967b453edSLaurent Vivier } else { 78067b453edSLaurent Vivier fp0 = floatx80_logn(a, status); 78167b453edSLaurent Vivier fp1 = packFloatx80(0, 0x3FFF, LIT64(0xB8AA3B295C17F0BC)); /* INV_L2 */ 78267b453edSLaurent Vivier 78367b453edSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 78467b453edSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 78567b453edSLaurent Vivier 78667b453edSLaurent Vivier a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L2 */ 78767b453edSLaurent Vivier } 78867b453edSLaurent Vivier 78967b453edSLaurent Vivier float_raise(float_flag_inexact, status); 79067b453edSLaurent Vivier 79167b453edSLaurent Vivier return a; 79267b453edSLaurent Vivier } 79340ad0873SLaurent Vivier 794*808d77bcSLucien Murray-Pitts /* 795*808d77bcSLucien Murray-Pitts * e to x 796*808d77bcSLucien Murray-Pitts */ 79740ad0873SLaurent Vivier 79840ad0873SLaurent Vivier floatx80 floatx80_etox(floatx80 a, float_status *status) 79940ad0873SLaurent Vivier { 80040ad0873SLaurent Vivier flag aSign; 80140ad0873SLaurent Vivier int32_t aExp; 80240ad0873SLaurent Vivier uint64_t aSig; 80340ad0873SLaurent Vivier 80440ad0873SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 80540ad0873SLaurent Vivier 80640ad0873SLaurent Vivier int32_t compact, n, j, k, m, m1; 80740ad0873SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, l2, scale, adjscale; 80840ad0873SLaurent Vivier flag adjflag; 80940ad0873SLaurent Vivier 81040ad0873SLaurent Vivier aSig = extractFloatx80Frac(a); 81140ad0873SLaurent Vivier aExp = extractFloatx80Exp(a); 81240ad0873SLaurent Vivier aSign = extractFloatx80Sign(a); 81340ad0873SLaurent Vivier 81440ad0873SLaurent Vivier if (aExp == 0x7FFF) { 81540ad0873SLaurent Vivier if ((uint64_t) (aSig << 1)) { 81640ad0873SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 81740ad0873SLaurent Vivier } 81840ad0873SLaurent Vivier if (aSign) { 81940ad0873SLaurent Vivier return packFloatx80(0, 0, 0); 82040ad0873SLaurent Vivier } 82140ad0873SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 82240ad0873SLaurent Vivier floatx80_infinity.low); 82340ad0873SLaurent Vivier } 82440ad0873SLaurent Vivier 82540ad0873SLaurent Vivier if (aExp == 0 && aSig == 0) { 82640ad0873SLaurent Vivier return packFloatx80(0, one_exp, one_sig); 82740ad0873SLaurent Vivier } 82840ad0873SLaurent Vivier 82940ad0873SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 83040ad0873SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 83140ad0873SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 83240ad0873SLaurent Vivier status->floatx80_rounding_precision = 80; 83340ad0873SLaurent Vivier 83440ad0873SLaurent Vivier adjflag = 0; 83540ad0873SLaurent Vivier 83640ad0873SLaurent Vivier if (aExp >= 0x3FBE) { /* |X| >= 2^(-65) */ 83740ad0873SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 83840ad0873SLaurent Vivier 83940ad0873SLaurent Vivier if (compact < 0x400CB167) { /* |X| < 16380 log2 */ 84040ad0873SLaurent Vivier fp0 = a; 84140ad0873SLaurent Vivier fp1 = a; 84240ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 84340ad0873SLaurent Vivier make_float32(0x42B8AA3B), status), 84440ad0873SLaurent Vivier status); /* 64/log2 * X */ 84540ad0873SLaurent Vivier adjflag = 0; 84640ad0873SLaurent Vivier n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */ 84740ad0873SLaurent Vivier fp0 = int32_to_floatx80(n, status); 84840ad0873SLaurent Vivier 84940ad0873SLaurent Vivier j = n & 0x3F; /* J = N mod 64 */ 85040ad0873SLaurent Vivier m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */ 85140ad0873SLaurent Vivier if (n < 0 && j) { 852*808d77bcSLucien Murray-Pitts /* 853*808d77bcSLucien Murray-Pitts * arithmetic right shift is division and 85440ad0873SLaurent Vivier * round towards minus infinity 85540ad0873SLaurent Vivier */ 85640ad0873SLaurent Vivier m--; 85740ad0873SLaurent Vivier } 85840ad0873SLaurent Vivier m += 0x3FFF; /* biased exponent of 2^(M) */ 85940ad0873SLaurent Vivier 86040ad0873SLaurent Vivier expcont1: 86140ad0873SLaurent Vivier fp2 = fp0; /* N */ 86240ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 86340ad0873SLaurent Vivier make_float32(0xBC317218), status), 86440ad0873SLaurent Vivier status); /* N * L1, L1 = lead(-log2/64) */ 86540ad0873SLaurent Vivier l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6)); 86640ad0873SLaurent Vivier fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */ 86740ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */ 86840ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* R */ 86940ad0873SLaurent Vivier 87040ad0873SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 87140ad0873SLaurent Vivier fp2 = float32_to_floatx80(make_float32(0x3AB60B70), 87240ad0873SLaurent Vivier status); /* A5 */ 87340ad0873SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A5 */ 87440ad0873SLaurent Vivier fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3C088895), 87540ad0873SLaurent Vivier status), fp1, 87640ad0873SLaurent Vivier status); /* fp3 is S*A4 */ 87740ad0873SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80(make_float64( 87840ad0873SLaurent Vivier 0x3FA5555555554431), status), 87940ad0873SLaurent Vivier status); /* fp2 is A3+S*A5 */ 88040ad0873SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80(make_float64( 88140ad0873SLaurent Vivier 0x3FC5555555554018), status), 88240ad0873SLaurent Vivier status); /* fp3 is A2+S*A4 */ 88340ad0873SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*(A3+S*A5) */ 88440ad0873SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* fp3 is S*(A2+S*A4) */ 88540ad0873SLaurent Vivier fp2 = floatx80_add(fp2, float32_to_floatx80( 88640ad0873SLaurent Vivier make_float32(0x3F000000), status), 88740ad0873SLaurent Vivier status); /* fp2 is A1+S*(A3+S*A5) */ 88840ad0873SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* fp3 IS R*S*(A2+S*A4) */ 88940ad0873SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, 89040ad0873SLaurent Vivier status); /* fp2 IS S*(A1+S*(A3+S*A5)) */ 89140ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp3, status); /* fp0 IS R+R*S*(A2+S*A4) */ 89240ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */ 89340ad0873SLaurent Vivier 89440ad0873SLaurent Vivier fp1 = exp_tbl[j]; 89540ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); /* 2^(J/64)*(Exp(R)-1) */ 89640ad0873SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status), 89740ad0873SLaurent Vivier status); /* accurate 2^(J/64) */ 89840ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp1, 89940ad0873SLaurent Vivier status); /* 2^(J/64) + 2^(J/64)*(Exp(R)-1) */ 90040ad0873SLaurent Vivier 90140ad0873SLaurent Vivier scale = packFloatx80(0, m, one_sig); 90240ad0873SLaurent Vivier if (adjflag) { 90340ad0873SLaurent Vivier adjscale = packFloatx80(0, m1, one_sig); 90440ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, adjscale, status); 90540ad0873SLaurent Vivier } 90640ad0873SLaurent Vivier 90740ad0873SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 90840ad0873SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 90940ad0873SLaurent Vivier 91040ad0873SLaurent Vivier a = floatx80_mul(fp0, scale, status); 91140ad0873SLaurent Vivier 91240ad0873SLaurent Vivier float_raise(float_flag_inexact, status); 91340ad0873SLaurent Vivier 91440ad0873SLaurent Vivier return a; 91540ad0873SLaurent Vivier } else { /* |X| >= 16380 log2 */ 91640ad0873SLaurent Vivier if (compact > 0x400CB27C) { /* |X| >= 16480 log2 */ 91740ad0873SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 91840ad0873SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 91940ad0873SLaurent Vivier if (aSign) { 92040ad0873SLaurent Vivier a = roundAndPackFloatx80( 92140ad0873SLaurent Vivier status->floatx80_rounding_precision, 92240ad0873SLaurent Vivier 0, -0x1000, aSig, 0, status); 92340ad0873SLaurent Vivier } else { 92440ad0873SLaurent Vivier a = roundAndPackFloatx80( 92540ad0873SLaurent Vivier status->floatx80_rounding_precision, 92640ad0873SLaurent Vivier 0, 0x8000, aSig, 0, status); 92740ad0873SLaurent Vivier } 92840ad0873SLaurent Vivier float_raise(float_flag_inexact, status); 92940ad0873SLaurent Vivier 93040ad0873SLaurent Vivier return a; 93140ad0873SLaurent Vivier } else { 93240ad0873SLaurent Vivier fp0 = a; 93340ad0873SLaurent Vivier fp1 = a; 93440ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 93540ad0873SLaurent Vivier make_float32(0x42B8AA3B), status), 93640ad0873SLaurent Vivier status); /* 64/log2 * X */ 93740ad0873SLaurent Vivier adjflag = 1; 93840ad0873SLaurent Vivier n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */ 93940ad0873SLaurent Vivier fp0 = int32_to_floatx80(n, status); 94040ad0873SLaurent Vivier 94140ad0873SLaurent Vivier j = n & 0x3F; /* J = N mod 64 */ 94240ad0873SLaurent Vivier /* NOTE: this is really arithmetic right shift by 6 */ 94340ad0873SLaurent Vivier k = n / 64; 94440ad0873SLaurent Vivier if (n < 0 && j) { 94540ad0873SLaurent Vivier /* arithmetic right shift is division and 94640ad0873SLaurent Vivier * round towards minus infinity 94740ad0873SLaurent Vivier */ 94840ad0873SLaurent Vivier k--; 94940ad0873SLaurent Vivier } 95040ad0873SLaurent Vivier /* NOTE: this is really arithmetic right shift by 1 */ 95140ad0873SLaurent Vivier m1 = k / 2; 95240ad0873SLaurent Vivier if (k < 0 && (k & 1)) { 95340ad0873SLaurent Vivier /* arithmetic right shift is division and 95440ad0873SLaurent Vivier * round towards minus infinity 95540ad0873SLaurent Vivier */ 95640ad0873SLaurent Vivier m1--; 95740ad0873SLaurent Vivier } 95840ad0873SLaurent Vivier m = k - m1; 95940ad0873SLaurent Vivier m1 += 0x3FFF; /* biased exponent of 2^(M1) */ 96040ad0873SLaurent Vivier m += 0x3FFF; /* biased exponent of 2^(M) */ 96140ad0873SLaurent Vivier 96240ad0873SLaurent Vivier goto expcont1; 96340ad0873SLaurent Vivier } 96440ad0873SLaurent Vivier } 96540ad0873SLaurent Vivier } else { /* |X| < 2^(-65) */ 96640ad0873SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 96740ad0873SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 96840ad0873SLaurent Vivier 96940ad0873SLaurent Vivier a = floatx80_add(a, float32_to_floatx80(make_float32(0x3F800000), 97040ad0873SLaurent Vivier status), status); /* 1 + X */ 97140ad0873SLaurent Vivier 97240ad0873SLaurent Vivier float_raise(float_flag_inexact, status); 97340ad0873SLaurent Vivier 97440ad0873SLaurent Vivier return a; 97540ad0873SLaurent Vivier } 97640ad0873SLaurent Vivier } 977068f1615SLaurent Vivier 978*808d77bcSLucien Murray-Pitts /* 979*808d77bcSLucien Murray-Pitts * 2 to x 980*808d77bcSLucien Murray-Pitts */ 981068f1615SLaurent Vivier 982068f1615SLaurent Vivier floatx80 floatx80_twotox(floatx80 a, float_status *status) 983068f1615SLaurent Vivier { 984068f1615SLaurent Vivier flag aSign; 985068f1615SLaurent Vivier int32_t aExp; 986068f1615SLaurent Vivier uint64_t aSig; 987068f1615SLaurent Vivier 988068f1615SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 989068f1615SLaurent Vivier 990068f1615SLaurent Vivier int32_t compact, n, j, l, m, m1; 991068f1615SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2; 992068f1615SLaurent Vivier 993068f1615SLaurent Vivier aSig = extractFloatx80Frac(a); 994068f1615SLaurent Vivier aExp = extractFloatx80Exp(a); 995068f1615SLaurent Vivier aSign = extractFloatx80Sign(a); 996068f1615SLaurent Vivier 997068f1615SLaurent Vivier if (aExp == 0x7FFF) { 998068f1615SLaurent Vivier if ((uint64_t) (aSig << 1)) { 999068f1615SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 1000068f1615SLaurent Vivier } 1001068f1615SLaurent Vivier if (aSign) { 1002068f1615SLaurent Vivier return packFloatx80(0, 0, 0); 1003068f1615SLaurent Vivier } 1004068f1615SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 1005068f1615SLaurent Vivier floatx80_infinity.low); 1006068f1615SLaurent Vivier } 1007068f1615SLaurent Vivier 1008068f1615SLaurent Vivier if (aExp == 0 && aSig == 0) { 1009068f1615SLaurent Vivier return packFloatx80(0, one_exp, one_sig); 1010068f1615SLaurent Vivier } 1011068f1615SLaurent Vivier 1012068f1615SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 1013068f1615SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 1014068f1615SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 1015068f1615SLaurent Vivier status->floatx80_rounding_precision = 80; 1016068f1615SLaurent Vivier 1017068f1615SLaurent Vivier fp0 = a; 1018068f1615SLaurent Vivier 1019068f1615SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 1020068f1615SLaurent Vivier 1021068f1615SLaurent Vivier if (compact < 0x3FB98000 || compact > 0x400D80C0) { 1022068f1615SLaurent Vivier /* |X| > 16480 or |X| < 2^(-70) */ 1023068f1615SLaurent Vivier if (compact > 0x3FFF8000) { /* |X| > 16480 */ 1024068f1615SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 1025068f1615SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 1026068f1615SLaurent Vivier 1027068f1615SLaurent Vivier if (aSign) { 1028068f1615SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 1029068f1615SLaurent Vivier 0, -0x1000, aSig, 0, status); 1030068f1615SLaurent Vivier } else { 1031068f1615SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 1032068f1615SLaurent Vivier 0, 0x8000, aSig, 0, status); 1033068f1615SLaurent Vivier } 1034068f1615SLaurent Vivier } else { /* |X| < 2^(-70) */ 1035068f1615SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 1036068f1615SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 1037068f1615SLaurent Vivier 1038068f1615SLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80( 1039068f1615SLaurent Vivier make_float32(0x3F800000), status), 1040068f1615SLaurent Vivier status); /* 1 + X */ 1041068f1615SLaurent Vivier 1042068f1615SLaurent Vivier float_raise(float_flag_inexact, status); 1043068f1615SLaurent Vivier 1044068f1615SLaurent Vivier return a; 1045068f1615SLaurent Vivier } 1046068f1615SLaurent Vivier } else { /* 2^(-70) <= |X| <= 16480 */ 1047068f1615SLaurent Vivier fp1 = fp0; /* X */ 1048068f1615SLaurent Vivier fp1 = floatx80_mul(fp1, float32_to_floatx80( 1049068f1615SLaurent Vivier make_float32(0x42800000), status), 1050068f1615SLaurent Vivier status); /* X * 64 */ 1051068f1615SLaurent Vivier n = floatx80_to_int32(fp1, status); 1052068f1615SLaurent Vivier fp1 = int32_to_floatx80(n, status); 1053068f1615SLaurent Vivier j = n & 0x3F; 1054068f1615SLaurent Vivier l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */ 1055068f1615SLaurent Vivier if (n < 0 && j) { 1056*808d77bcSLucien Murray-Pitts /* 1057*808d77bcSLucien Murray-Pitts * arithmetic right shift is division and 1058068f1615SLaurent Vivier * round towards minus infinity 1059068f1615SLaurent Vivier */ 1060068f1615SLaurent Vivier l--; 1061068f1615SLaurent Vivier } 1062068f1615SLaurent Vivier m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */ 1063068f1615SLaurent Vivier if (l < 0 && (l & 1)) { 1064*808d77bcSLucien Murray-Pitts /* 1065*808d77bcSLucien Murray-Pitts * arithmetic right shift is division and 1066068f1615SLaurent Vivier * round towards minus infinity 1067068f1615SLaurent Vivier */ 1068068f1615SLaurent Vivier m--; 1069068f1615SLaurent Vivier } 1070068f1615SLaurent Vivier m1 = l - m; 1071068f1615SLaurent Vivier m1 += 0x3FFF; /* ADJFACT IS 2^(M') */ 1072068f1615SLaurent Vivier 1073068f1615SLaurent Vivier adjfact = packFloatx80(0, m1, one_sig); 1074068f1615SLaurent Vivier fact1 = exp2_tbl[j]; 1075068f1615SLaurent Vivier fact1.high += m; 1076068f1615SLaurent Vivier fact2.high = exp2_tbl2[j] >> 16; 1077068f1615SLaurent Vivier fact2.high += m; 1078068f1615SLaurent Vivier fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF); 1079068f1615SLaurent Vivier fact2.low <<= 48; 1080068f1615SLaurent Vivier 1081068f1615SLaurent Vivier fp1 = floatx80_mul(fp1, float32_to_floatx80( 1082068f1615SLaurent Vivier make_float32(0x3C800000), status), 1083068f1615SLaurent Vivier status); /* (1/64)*N */ 1084068f1615SLaurent Vivier fp0 = floatx80_sub(fp0, fp1, status); /* X - (1/64)*INT(64 X) */ 1085068f1615SLaurent Vivier fp2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); /* LOG2 */ 1086068f1615SLaurent Vivier fp0 = floatx80_mul(fp0, fp2, status); /* R */ 1087068f1615SLaurent Vivier 1088068f1615SLaurent Vivier /* EXPR */ 1089068f1615SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 1090068f1615SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2), 1091068f1615SLaurent Vivier status); /* A5 */ 1092068f1615SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C), 1093068f1615SLaurent Vivier status); /* A4 */ 1094068f1615SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */ 1095068f1615SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */ 1096068f1615SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 1097068f1615SLaurent Vivier make_float64(0x3FA5555555554CC1), status), 1098068f1615SLaurent Vivier status); /* A3+S*A5 */ 1099068f1615SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 1100068f1615SLaurent Vivier make_float64(0x3FC5555555554A54), status), 1101068f1615SLaurent Vivier status); /* A2+S*A4 */ 1102068f1615SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */ 1103068f1615SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */ 1104068f1615SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 1105068f1615SLaurent Vivier make_float64(0x3FE0000000000000), status), 1106068f1615SLaurent Vivier status); /* A1+S*(A3+S*A5) */ 1107068f1615SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */ 1108068f1615SLaurent Vivier 1109068f1615SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */ 1110068f1615SLaurent Vivier fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */ 1111068f1615SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */ 1112068f1615SLaurent Vivier 1113068f1615SLaurent Vivier fp0 = floatx80_mul(fp0, fact1, status); 1114068f1615SLaurent Vivier fp0 = floatx80_add(fp0, fact2, status); 1115068f1615SLaurent Vivier fp0 = floatx80_add(fp0, fact1, status); 1116068f1615SLaurent Vivier 1117068f1615SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 1118068f1615SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 1119068f1615SLaurent Vivier 1120068f1615SLaurent Vivier a = floatx80_mul(fp0, adjfact, status); 1121068f1615SLaurent Vivier 1122068f1615SLaurent Vivier float_raise(float_flag_inexact, status); 1123068f1615SLaurent Vivier 1124068f1615SLaurent Vivier return a; 1125068f1615SLaurent Vivier } 1126068f1615SLaurent Vivier } 11276c25be6eSLaurent Vivier 1128*808d77bcSLucien Murray-Pitts /* 1129*808d77bcSLucien Murray-Pitts * 10 to x 1130*808d77bcSLucien Murray-Pitts */ 11316c25be6eSLaurent Vivier 11326c25be6eSLaurent Vivier floatx80 floatx80_tentox(floatx80 a, float_status *status) 11336c25be6eSLaurent Vivier { 11346c25be6eSLaurent Vivier flag aSign; 11356c25be6eSLaurent Vivier int32_t aExp; 11366c25be6eSLaurent Vivier uint64_t aSig; 11376c25be6eSLaurent Vivier 11386c25be6eSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 11396c25be6eSLaurent Vivier 11406c25be6eSLaurent Vivier int32_t compact, n, j, l, m, m1; 11416c25be6eSLaurent Vivier floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2; 11426c25be6eSLaurent Vivier 11436c25be6eSLaurent Vivier aSig = extractFloatx80Frac(a); 11446c25be6eSLaurent Vivier aExp = extractFloatx80Exp(a); 11456c25be6eSLaurent Vivier aSign = extractFloatx80Sign(a); 11466c25be6eSLaurent Vivier 11476c25be6eSLaurent Vivier if (aExp == 0x7FFF) { 11486c25be6eSLaurent Vivier if ((uint64_t) (aSig << 1)) { 11496c25be6eSLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 11506c25be6eSLaurent Vivier } 11516c25be6eSLaurent Vivier if (aSign) { 11526c25be6eSLaurent Vivier return packFloatx80(0, 0, 0); 11536c25be6eSLaurent Vivier } 11546c25be6eSLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 11556c25be6eSLaurent Vivier floatx80_infinity.low); 11566c25be6eSLaurent Vivier } 11576c25be6eSLaurent Vivier 11586c25be6eSLaurent Vivier if (aExp == 0 && aSig == 0) { 11596c25be6eSLaurent Vivier return packFloatx80(0, one_exp, one_sig); 11606c25be6eSLaurent Vivier } 11616c25be6eSLaurent Vivier 11626c25be6eSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 11636c25be6eSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 11646c25be6eSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 11656c25be6eSLaurent Vivier status->floatx80_rounding_precision = 80; 11666c25be6eSLaurent Vivier 11676c25be6eSLaurent Vivier fp0 = a; 11686c25be6eSLaurent Vivier 11696c25be6eSLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 11706c25be6eSLaurent Vivier 11716c25be6eSLaurent Vivier if (compact < 0x3FB98000 || compact > 0x400B9B07) { 11726c25be6eSLaurent Vivier /* |X| > 16480 LOG2/LOG10 or |X| < 2^(-70) */ 11736c25be6eSLaurent Vivier if (compact > 0x3FFF8000) { /* |X| > 16480 */ 11746c25be6eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 11756c25be6eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 11766c25be6eSLaurent Vivier 11776c25be6eSLaurent Vivier if (aSign) { 11786c25be6eSLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 11796c25be6eSLaurent Vivier 0, -0x1000, aSig, 0, status); 11806c25be6eSLaurent Vivier } else { 11816c25be6eSLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 11826c25be6eSLaurent Vivier 0, 0x8000, aSig, 0, status); 11836c25be6eSLaurent Vivier } 11846c25be6eSLaurent Vivier } else { /* |X| < 2^(-70) */ 11856c25be6eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 11866c25be6eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 11876c25be6eSLaurent Vivier 11886c25be6eSLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80( 11896c25be6eSLaurent Vivier make_float32(0x3F800000), status), 11906c25be6eSLaurent Vivier status); /* 1 + X */ 11916c25be6eSLaurent Vivier 11926c25be6eSLaurent Vivier float_raise(float_flag_inexact, status); 11936c25be6eSLaurent Vivier 11946c25be6eSLaurent Vivier return a; 11956c25be6eSLaurent Vivier } 11966c25be6eSLaurent Vivier } else { /* 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10 */ 11976c25be6eSLaurent Vivier fp1 = fp0; /* X */ 11986c25be6eSLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 11996c25be6eSLaurent Vivier make_float64(0x406A934F0979A371), 12006c25be6eSLaurent Vivier status), status); /* X*64*LOG10/LOG2 */ 12016c25be6eSLaurent Vivier n = floatx80_to_int32(fp1, status); /* N=INT(X*64*LOG10/LOG2) */ 12026c25be6eSLaurent Vivier fp1 = int32_to_floatx80(n, status); 12036c25be6eSLaurent Vivier 12046c25be6eSLaurent Vivier j = n & 0x3F; 12056c25be6eSLaurent Vivier l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */ 12066c25be6eSLaurent Vivier if (n < 0 && j) { 1207*808d77bcSLucien Murray-Pitts /* 1208*808d77bcSLucien Murray-Pitts * arithmetic right shift is division and 12096c25be6eSLaurent Vivier * round towards minus infinity 12106c25be6eSLaurent Vivier */ 12116c25be6eSLaurent Vivier l--; 12126c25be6eSLaurent Vivier } 12136c25be6eSLaurent Vivier m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */ 12146c25be6eSLaurent Vivier if (l < 0 && (l & 1)) { 1215*808d77bcSLucien Murray-Pitts /* 1216*808d77bcSLucien Murray-Pitts * arithmetic right shift is division and 12176c25be6eSLaurent Vivier * round towards minus infinity 12186c25be6eSLaurent Vivier */ 12196c25be6eSLaurent Vivier m--; 12206c25be6eSLaurent Vivier } 12216c25be6eSLaurent Vivier m1 = l - m; 12226c25be6eSLaurent Vivier m1 += 0x3FFF; /* ADJFACT IS 2^(M') */ 12236c25be6eSLaurent Vivier 12246c25be6eSLaurent Vivier adjfact = packFloatx80(0, m1, one_sig); 12256c25be6eSLaurent Vivier fact1 = exp2_tbl[j]; 12266c25be6eSLaurent Vivier fact1.high += m; 12276c25be6eSLaurent Vivier fact2.high = exp2_tbl2[j] >> 16; 12286c25be6eSLaurent Vivier fact2.high += m; 12296c25be6eSLaurent Vivier fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF); 12306c25be6eSLaurent Vivier fact2.low <<= 48; 12316c25be6eSLaurent Vivier 12326c25be6eSLaurent Vivier fp2 = fp1; /* N */ 12336c25be6eSLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 12346c25be6eSLaurent Vivier make_float64(0x3F734413509F8000), status), 12356c25be6eSLaurent Vivier status); /* N*(LOG2/64LOG10)_LEAD */ 12366c25be6eSLaurent Vivier fp3 = packFloatx80(1, 0x3FCD, LIT64(0xC0219DC1DA994FD2)); 12376c25be6eSLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* N*(LOG2/64LOG10)_TRAIL */ 12386c25be6eSLaurent Vivier fp0 = floatx80_sub(fp0, fp1, status); /* X - N L_LEAD */ 12396c25be6eSLaurent Vivier fp0 = floatx80_sub(fp0, fp2, status); /* X - N L_TRAIL */ 12406c25be6eSLaurent Vivier fp2 = packFloatx80(0, 0x4000, LIT64(0x935D8DDDAAA8AC17)); /* LOG10 */ 12416c25be6eSLaurent Vivier fp0 = floatx80_mul(fp0, fp2, status); /* R */ 12426c25be6eSLaurent Vivier 12436c25be6eSLaurent Vivier /* EXPR */ 12446c25be6eSLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 12456c25be6eSLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2), 12466c25be6eSLaurent Vivier status); /* A5 */ 12476c25be6eSLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C), 12486c25be6eSLaurent Vivier status); /* A4 */ 12496c25be6eSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */ 12506c25be6eSLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */ 12516c25be6eSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 12526c25be6eSLaurent Vivier make_float64(0x3FA5555555554CC1), status), 12536c25be6eSLaurent Vivier status); /* A3+S*A5 */ 12546c25be6eSLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 12556c25be6eSLaurent Vivier make_float64(0x3FC5555555554A54), status), 12566c25be6eSLaurent Vivier status); /* A2+S*A4 */ 12576c25be6eSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */ 12586c25be6eSLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */ 12596c25be6eSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 12606c25be6eSLaurent Vivier make_float64(0x3FE0000000000000), status), 12616c25be6eSLaurent Vivier status); /* A1+S*(A3+S*A5) */ 12626c25be6eSLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */ 12636c25be6eSLaurent Vivier 12646c25be6eSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */ 12656c25be6eSLaurent Vivier fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */ 12666c25be6eSLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */ 12676c25be6eSLaurent Vivier 12686c25be6eSLaurent Vivier fp0 = floatx80_mul(fp0, fact1, status); 12696c25be6eSLaurent Vivier fp0 = floatx80_add(fp0, fact2, status); 12706c25be6eSLaurent Vivier fp0 = floatx80_add(fp0, fact1, status); 12716c25be6eSLaurent Vivier 12726c25be6eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 12736c25be6eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 12746c25be6eSLaurent Vivier 12756c25be6eSLaurent Vivier a = floatx80_mul(fp0, adjfact, status); 12766c25be6eSLaurent Vivier 12776c25be6eSLaurent Vivier float_raise(float_flag_inexact, status); 12786c25be6eSLaurent Vivier 12796c25be6eSLaurent Vivier return a; 12806c25be6eSLaurent Vivier } 12816c25be6eSLaurent Vivier } 128227340180SLaurent Vivier 1283*808d77bcSLucien Murray-Pitts /* 1284*808d77bcSLucien Murray-Pitts * Tangent 1285*808d77bcSLucien Murray-Pitts */ 128627340180SLaurent Vivier 128727340180SLaurent Vivier floatx80 floatx80_tan(floatx80 a, float_status *status) 128827340180SLaurent Vivier { 128927340180SLaurent Vivier flag aSign, xSign; 129027340180SLaurent Vivier int32_t aExp, xExp; 129127340180SLaurent Vivier uint64_t aSig, xSig; 129227340180SLaurent Vivier 129327340180SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 129427340180SLaurent Vivier 129527340180SLaurent Vivier int32_t compact, l, n, j; 129627340180SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, fp4, fp5, invtwopi, twopi1, twopi2; 129727340180SLaurent Vivier float32 twoto63; 129827340180SLaurent Vivier flag endflag; 129927340180SLaurent Vivier 130027340180SLaurent Vivier aSig = extractFloatx80Frac(a); 130127340180SLaurent Vivier aExp = extractFloatx80Exp(a); 130227340180SLaurent Vivier aSign = extractFloatx80Sign(a); 130327340180SLaurent Vivier 130427340180SLaurent Vivier if (aExp == 0x7FFF) { 130527340180SLaurent Vivier if ((uint64_t) (aSig << 1)) { 130627340180SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 130727340180SLaurent Vivier } 130827340180SLaurent Vivier float_raise(float_flag_invalid, status); 130927340180SLaurent Vivier return floatx80_default_nan(status); 131027340180SLaurent Vivier } 131127340180SLaurent Vivier 131227340180SLaurent Vivier if (aExp == 0 && aSig == 0) { 131327340180SLaurent Vivier return packFloatx80(aSign, 0, 0); 131427340180SLaurent Vivier } 131527340180SLaurent Vivier 131627340180SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 131727340180SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 131827340180SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 131927340180SLaurent Vivier status->floatx80_rounding_precision = 80; 132027340180SLaurent Vivier 132127340180SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 132227340180SLaurent Vivier 132327340180SLaurent Vivier fp0 = a; 132427340180SLaurent Vivier 132527340180SLaurent Vivier if (compact < 0x3FD78000 || compact > 0x4004BC7E) { 132627340180SLaurent Vivier /* 2^(-40) > |X| > 15 PI */ 132727340180SLaurent Vivier if (compact > 0x3FFF8000) { /* |X| >= 15 PI */ 132827340180SLaurent Vivier /* REDUCEX */ 132927340180SLaurent Vivier fp1 = packFloatx80(0, 0, 0); 133027340180SLaurent Vivier if (compact == 0x7FFEFFFF) { 133127340180SLaurent Vivier twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, 133227340180SLaurent Vivier LIT64(0xC90FDAA200000000)); 133327340180SLaurent Vivier twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, 133427340180SLaurent Vivier LIT64(0x85A308D300000000)); 133527340180SLaurent Vivier fp0 = floatx80_add(fp0, twopi1, status); 133627340180SLaurent Vivier fp1 = fp0; 133727340180SLaurent Vivier fp0 = floatx80_add(fp0, twopi2, status); 133827340180SLaurent Vivier fp1 = floatx80_sub(fp1, fp0, status); 133927340180SLaurent Vivier fp1 = floatx80_add(fp1, twopi2, status); 134027340180SLaurent Vivier } 134127340180SLaurent Vivier loop: 134227340180SLaurent Vivier xSign = extractFloatx80Sign(fp0); 134327340180SLaurent Vivier xExp = extractFloatx80Exp(fp0); 134427340180SLaurent Vivier xExp -= 0x3FFF; 134527340180SLaurent Vivier if (xExp <= 28) { 134627340180SLaurent Vivier l = 0; 134727340180SLaurent Vivier endflag = 1; 134827340180SLaurent Vivier } else { 134927340180SLaurent Vivier l = xExp - 27; 135027340180SLaurent Vivier endflag = 0; 135127340180SLaurent Vivier } 135227340180SLaurent Vivier invtwopi = packFloatx80(0, 0x3FFE - l, 135327340180SLaurent Vivier LIT64(0xA2F9836E4E44152A)); /* INVTWOPI */ 135427340180SLaurent Vivier twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); 135527340180SLaurent Vivier twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); 135627340180SLaurent Vivier 135727340180SLaurent Vivier /* SIGN(INARG)*2^63 IN SGL */ 135827340180SLaurent Vivier twoto63 = packFloat32(xSign, 0xBE, 0); 135927340180SLaurent Vivier 136027340180SLaurent Vivier fp2 = floatx80_mul(fp0, invtwopi, status); 136127340180SLaurent Vivier fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), 136227340180SLaurent Vivier status); /* THE FRACT PART OF FP2 IS ROUNDED */ 136327340180SLaurent Vivier fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), 136427340180SLaurent Vivier status); /* FP2 is N */ 136527340180SLaurent Vivier fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */ 136627340180SLaurent Vivier fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */ 136727340180SLaurent Vivier fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */ 136827340180SLaurent Vivier fp4 = floatx80_sub(fp4, fp3, status); /* W-P */ 136927340180SLaurent Vivier fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */ 137027340180SLaurent Vivier fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */ 137127340180SLaurent Vivier fp3 = fp0; /* FP3 is A */ 137227340180SLaurent Vivier fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */ 137327340180SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */ 137427340180SLaurent Vivier 137527340180SLaurent Vivier if (endflag > 0) { 137627340180SLaurent Vivier n = floatx80_to_int32(fp2, status); 137727340180SLaurent Vivier goto tancont; 137827340180SLaurent Vivier } 137927340180SLaurent Vivier fp3 = floatx80_sub(fp3, fp0, status); /* A-R */ 138027340180SLaurent Vivier fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */ 138127340180SLaurent Vivier goto loop; 138227340180SLaurent Vivier } else { 138327340180SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 138427340180SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 138527340180SLaurent Vivier 138627340180SLaurent Vivier a = floatx80_move(a, status); 138727340180SLaurent Vivier 138827340180SLaurent Vivier float_raise(float_flag_inexact, status); 138927340180SLaurent Vivier 139027340180SLaurent Vivier return a; 139127340180SLaurent Vivier } 139227340180SLaurent Vivier } else { 139327340180SLaurent Vivier fp1 = floatx80_mul(fp0, float64_to_floatx80( 139427340180SLaurent Vivier make_float64(0x3FE45F306DC9C883), status), 139527340180SLaurent Vivier status); /* X*2/PI */ 139627340180SLaurent Vivier 139727340180SLaurent Vivier n = floatx80_to_int32(fp1, status); 139827340180SLaurent Vivier j = 32 + n; 139927340180SLaurent Vivier 140027340180SLaurent Vivier fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */ 140127340180SLaurent Vivier fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status), 140227340180SLaurent Vivier status); /* FP0 IS R = (X-Y1)-Y2 */ 140327340180SLaurent Vivier 140427340180SLaurent Vivier tancont: 140527340180SLaurent Vivier if (n & 1) { 140627340180SLaurent Vivier /* NODD */ 140727340180SLaurent Vivier fp1 = fp0; /* R */ 140827340180SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 140927340180SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688), 141027340180SLaurent Vivier status); /* Q4 */ 141127340180SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04), 141227340180SLaurent Vivier status); /* P3 */ 141327340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* SQ4 */ 141427340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* SP3 */ 141527340180SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 141627340180SLaurent Vivier make_float64(0xBF346F59B39BA65F), status), 141727340180SLaurent Vivier status); /* Q3+SQ4 */ 141827340180SLaurent Vivier fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00)); 141927340180SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */ 142027340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* S(Q3+SQ4) */ 142127340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* S(P2+SP3) */ 142227340180SLaurent Vivier fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1)); 142327340180SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */ 142427340180SLaurent Vivier fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA)); 142527340180SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */ 142627340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* S(Q2+S(Q3+SQ4)) */ 142727340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* S(P1+S(P2+SP3)) */ 142827340180SLaurent Vivier fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE)); 142927340180SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */ 143027340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* RS(P1+S(P2+SP3)) */ 143127340180SLaurent Vivier fp0 = floatx80_mul(fp0, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */ 143227340180SLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); /* R+RS(P1+S(P2+SP3)) */ 143327340180SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80( 143427340180SLaurent Vivier make_float32(0x3F800000), status), 143527340180SLaurent Vivier status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */ 143627340180SLaurent Vivier 143727340180SLaurent Vivier xSign = extractFloatx80Sign(fp1); 143827340180SLaurent Vivier xExp = extractFloatx80Exp(fp1); 143927340180SLaurent Vivier xSig = extractFloatx80Frac(fp1); 144027340180SLaurent Vivier xSign ^= 1; 144127340180SLaurent Vivier fp1 = packFloatx80(xSign, xExp, xSig); 144227340180SLaurent Vivier 144327340180SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 144427340180SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 144527340180SLaurent Vivier 144627340180SLaurent Vivier a = floatx80_div(fp0, fp1, status); 144727340180SLaurent Vivier 144827340180SLaurent Vivier float_raise(float_flag_inexact, status); 144927340180SLaurent Vivier 145027340180SLaurent Vivier return a; 145127340180SLaurent Vivier } else { 145227340180SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 145327340180SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688), 145427340180SLaurent Vivier status); /* Q4 */ 145527340180SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04), 145627340180SLaurent Vivier status); /* P3 */ 145727340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* SQ4 */ 145827340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* SP3 */ 145927340180SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 146027340180SLaurent Vivier make_float64(0xBF346F59B39BA65F), status), 146127340180SLaurent Vivier status); /* Q3+SQ4 */ 146227340180SLaurent Vivier fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00)); 146327340180SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */ 146427340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S(Q3+SQ4) */ 146527340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S(P2+SP3) */ 146627340180SLaurent Vivier fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1)); 146727340180SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */ 146827340180SLaurent Vivier fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA)); 146927340180SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */ 147027340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S(Q2+S(Q3+SQ4)) */ 147127340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S(P1+S(P2+SP3)) */ 147227340180SLaurent Vivier fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE)); 147327340180SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */ 147427340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* RS(P1+S(P2+SP3)) */ 147527340180SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */ 147627340180SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* R+RS(P1+S(P2+SP3)) */ 147727340180SLaurent Vivier fp1 = floatx80_add(fp1, float32_to_floatx80( 147827340180SLaurent Vivier make_float32(0x3F800000), status), 147927340180SLaurent Vivier status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */ 148027340180SLaurent Vivier 148127340180SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 148227340180SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 148327340180SLaurent Vivier 148427340180SLaurent Vivier a = floatx80_div(fp0, fp1, status); 148527340180SLaurent Vivier 148627340180SLaurent Vivier float_raise(float_flag_inexact, status); 148727340180SLaurent Vivier 148827340180SLaurent Vivier return a; 148927340180SLaurent Vivier } 149027340180SLaurent Vivier } 149127340180SLaurent Vivier } 14925add1ac4SLaurent Vivier 1493*808d77bcSLucien Murray-Pitts /* 1494*808d77bcSLucien Murray-Pitts * Sine 1495*808d77bcSLucien Murray-Pitts */ 14965add1ac4SLaurent Vivier 14975add1ac4SLaurent Vivier floatx80 floatx80_sin(floatx80 a, float_status *status) 14985add1ac4SLaurent Vivier { 14995add1ac4SLaurent Vivier flag aSign, xSign; 15005add1ac4SLaurent Vivier int32_t aExp, xExp; 15015add1ac4SLaurent Vivier uint64_t aSig, xSig; 15025add1ac4SLaurent Vivier 15035add1ac4SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 15045add1ac4SLaurent Vivier 15055add1ac4SLaurent Vivier int32_t compact, l, n, j; 15065add1ac4SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2; 15075add1ac4SLaurent Vivier float32 posneg1, twoto63; 15086361d298SLaurent Vivier flag endflag; 15095add1ac4SLaurent Vivier 15105add1ac4SLaurent Vivier aSig = extractFloatx80Frac(a); 15115add1ac4SLaurent Vivier aExp = extractFloatx80Exp(a); 15125add1ac4SLaurent Vivier aSign = extractFloatx80Sign(a); 15135add1ac4SLaurent Vivier 15145add1ac4SLaurent Vivier if (aExp == 0x7FFF) { 15155add1ac4SLaurent Vivier if ((uint64_t) (aSig << 1)) { 15165add1ac4SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 15175add1ac4SLaurent Vivier } 15185add1ac4SLaurent Vivier float_raise(float_flag_invalid, status); 15195add1ac4SLaurent Vivier return floatx80_default_nan(status); 15205add1ac4SLaurent Vivier } 15215add1ac4SLaurent Vivier 15225add1ac4SLaurent Vivier if (aExp == 0 && aSig == 0) { 15235add1ac4SLaurent Vivier return packFloatx80(aSign, 0, 0); 15245add1ac4SLaurent Vivier } 15255add1ac4SLaurent Vivier 15265add1ac4SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 15275add1ac4SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 15285add1ac4SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 15295add1ac4SLaurent Vivier status->floatx80_rounding_precision = 80; 15305add1ac4SLaurent Vivier 15315add1ac4SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 15325add1ac4SLaurent Vivier 15335add1ac4SLaurent Vivier fp0 = a; 15345add1ac4SLaurent Vivier 15355add1ac4SLaurent Vivier if (compact < 0x3FD78000 || compact > 0x4004BC7E) { 15365add1ac4SLaurent Vivier /* 2^(-40) > |X| > 15 PI */ 15375add1ac4SLaurent Vivier if (compact > 0x3FFF8000) { /* |X| >= 15 PI */ 15385add1ac4SLaurent Vivier /* REDUCEX */ 15395add1ac4SLaurent Vivier fp1 = packFloatx80(0, 0, 0); 15405add1ac4SLaurent Vivier if (compact == 0x7FFEFFFF) { 15415add1ac4SLaurent Vivier twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, 15425add1ac4SLaurent Vivier LIT64(0xC90FDAA200000000)); 15435add1ac4SLaurent Vivier twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, 15445add1ac4SLaurent Vivier LIT64(0x85A308D300000000)); 15455add1ac4SLaurent Vivier fp0 = floatx80_add(fp0, twopi1, status); 15465add1ac4SLaurent Vivier fp1 = fp0; 15475add1ac4SLaurent Vivier fp0 = floatx80_add(fp0, twopi2, status); 15485add1ac4SLaurent Vivier fp1 = floatx80_sub(fp1, fp0, status); 15495add1ac4SLaurent Vivier fp1 = floatx80_add(fp1, twopi2, status); 15505add1ac4SLaurent Vivier } 15515add1ac4SLaurent Vivier loop: 15525add1ac4SLaurent Vivier xSign = extractFloatx80Sign(fp0); 15535add1ac4SLaurent Vivier xExp = extractFloatx80Exp(fp0); 15545add1ac4SLaurent Vivier xExp -= 0x3FFF; 15555add1ac4SLaurent Vivier if (xExp <= 28) { 15565add1ac4SLaurent Vivier l = 0; 15575add1ac4SLaurent Vivier endflag = 1; 15585add1ac4SLaurent Vivier } else { 15595add1ac4SLaurent Vivier l = xExp - 27; 15605add1ac4SLaurent Vivier endflag = 0; 15615add1ac4SLaurent Vivier } 15625add1ac4SLaurent Vivier invtwopi = packFloatx80(0, 0x3FFE - l, 15635add1ac4SLaurent Vivier LIT64(0xA2F9836E4E44152A)); /* INVTWOPI */ 15645add1ac4SLaurent Vivier twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); 15655add1ac4SLaurent Vivier twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); 15665add1ac4SLaurent Vivier 15675add1ac4SLaurent Vivier /* SIGN(INARG)*2^63 IN SGL */ 15685add1ac4SLaurent Vivier twoto63 = packFloat32(xSign, 0xBE, 0); 15695add1ac4SLaurent Vivier 15705add1ac4SLaurent Vivier fp2 = floatx80_mul(fp0, invtwopi, status); 15715add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), 15725add1ac4SLaurent Vivier status); /* THE FRACT PART OF FP2 IS ROUNDED */ 15735add1ac4SLaurent Vivier fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), 15745add1ac4SLaurent Vivier status); /* FP2 is N */ 15755add1ac4SLaurent Vivier fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */ 15765add1ac4SLaurent Vivier fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */ 15775add1ac4SLaurent Vivier fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */ 15785add1ac4SLaurent Vivier fp4 = floatx80_sub(fp4, fp3, status); /* W-P */ 15795add1ac4SLaurent Vivier fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */ 15805add1ac4SLaurent Vivier fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */ 15815add1ac4SLaurent Vivier fp3 = fp0; /* FP3 is A */ 15825add1ac4SLaurent Vivier fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */ 15835add1ac4SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */ 15845add1ac4SLaurent Vivier 15855add1ac4SLaurent Vivier if (endflag > 0) { 15865add1ac4SLaurent Vivier n = floatx80_to_int32(fp2, status); 15875add1ac4SLaurent Vivier goto sincont; 15885add1ac4SLaurent Vivier } 15895add1ac4SLaurent Vivier fp3 = floatx80_sub(fp3, fp0, status); /* A-R */ 15905add1ac4SLaurent Vivier fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */ 15915add1ac4SLaurent Vivier goto loop; 15925add1ac4SLaurent Vivier } else { 15935add1ac4SLaurent Vivier /* SINSM */ 15945add1ac4SLaurent Vivier fp0 = float32_to_floatx80(make_float32(0x3F800000), 15955add1ac4SLaurent Vivier status); /* 1 */ 15965add1ac4SLaurent Vivier 15975add1ac4SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 15985add1ac4SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 15995add1ac4SLaurent Vivier 16005add1ac4SLaurent Vivier /* SINTINY */ 16015add1ac4SLaurent Vivier a = floatx80_move(a, status); 16025add1ac4SLaurent Vivier float_raise(float_flag_inexact, status); 16035add1ac4SLaurent Vivier 16045add1ac4SLaurent Vivier return a; 16055add1ac4SLaurent Vivier } 16065add1ac4SLaurent Vivier } else { 16075add1ac4SLaurent Vivier fp1 = floatx80_mul(fp0, float64_to_floatx80( 16085add1ac4SLaurent Vivier make_float64(0x3FE45F306DC9C883), status), 16095add1ac4SLaurent Vivier status); /* X*2/PI */ 16105add1ac4SLaurent Vivier 16115add1ac4SLaurent Vivier n = floatx80_to_int32(fp1, status); 16125add1ac4SLaurent Vivier j = 32 + n; 16135add1ac4SLaurent Vivier 16145add1ac4SLaurent Vivier fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */ 16155add1ac4SLaurent Vivier fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status), 16165add1ac4SLaurent Vivier status); /* FP0 IS R = (X-Y1)-Y2 */ 16175add1ac4SLaurent Vivier 16185add1ac4SLaurent Vivier sincont: 16196361d298SLaurent Vivier if (n & 1) { 16205add1ac4SLaurent Vivier /* COSPOLY */ 16215add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */ 16225add1ac4SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */ 16235add1ac4SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3), 16245add1ac4SLaurent Vivier status); /* B8 */ 16255add1ac4SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19), 16265add1ac4SLaurent Vivier status); /* B7 */ 16275add1ac4SLaurent Vivier 16285add1ac4SLaurent Vivier xSign = extractFloatx80Sign(fp0); /* X IS S */ 16295add1ac4SLaurent Vivier xExp = extractFloatx80Exp(fp0); 16305add1ac4SLaurent Vivier xSig = extractFloatx80Frac(fp0); 16315add1ac4SLaurent Vivier 16326361d298SLaurent Vivier if ((n >> 1) & 1) { 16335add1ac4SLaurent Vivier xSign ^= 1; 16345add1ac4SLaurent Vivier posneg1 = make_float32(0xBF800000); /* -1 */ 16355add1ac4SLaurent Vivier } else { 16365add1ac4SLaurent Vivier xSign ^= 0; 16375add1ac4SLaurent Vivier posneg1 = make_float32(0x3F800000); /* 1 */ 16385add1ac4SLaurent Vivier } /* X IS NOW R'= SGN*R */ 16395add1ac4SLaurent Vivier 16405add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */ 16415add1ac4SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */ 16425add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 16435add1ac4SLaurent Vivier make_float64(0x3E21EED90612C972), status), 16445add1ac4SLaurent Vivier status); /* B6+TB8 */ 16455add1ac4SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 16465add1ac4SLaurent Vivier make_float64(0xBE927E4FB79D9FCF), status), 16475add1ac4SLaurent Vivier status); /* B5+TB7 */ 16485add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */ 16495add1ac4SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */ 16505add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 16515add1ac4SLaurent Vivier make_float64(0x3EFA01A01A01D423), status), 16525add1ac4SLaurent Vivier status); /* B4+T(B6+TB8) */ 16535add1ac4SLaurent Vivier fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438)); 16545add1ac4SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */ 16555add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */ 16565add1ac4SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */ 16575add1ac4SLaurent Vivier fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E)); 16585add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */ 16595add1ac4SLaurent Vivier fp1 = floatx80_add(fp1, float32_to_floatx80( 16605add1ac4SLaurent Vivier make_float32(0xBF000000), status), 16615add1ac4SLaurent Vivier status); /* B1+T(B3+T(B5+TB7)) */ 16625add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */ 16635add1ac4SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* [B1+T(B3+T(B5+TB7))]+ 16645add1ac4SLaurent Vivier * [S(B2+T(B4+T(B6+TB8)))] 16655add1ac4SLaurent Vivier */ 16665add1ac4SLaurent Vivier 16675add1ac4SLaurent Vivier x = packFloatx80(xSign, xExp, xSig); 16685add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, x, status); 16695add1ac4SLaurent Vivier 16705add1ac4SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 16715add1ac4SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 16725add1ac4SLaurent Vivier 16735add1ac4SLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status); 16745add1ac4SLaurent Vivier 16755add1ac4SLaurent Vivier float_raise(float_flag_inexact, status); 16765add1ac4SLaurent Vivier 16775add1ac4SLaurent Vivier return a; 16785add1ac4SLaurent Vivier } else { 16795add1ac4SLaurent Vivier /* SINPOLY */ 16805add1ac4SLaurent Vivier xSign = extractFloatx80Sign(fp0); /* X IS R */ 16815add1ac4SLaurent Vivier xExp = extractFloatx80Exp(fp0); 16825add1ac4SLaurent Vivier xSig = extractFloatx80Frac(fp0); 16835add1ac4SLaurent Vivier 16846361d298SLaurent Vivier xSign ^= (n >> 1) & 1; /* X IS NOW R'= SGN*R */ 16855add1ac4SLaurent Vivier 16865add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */ 16875add1ac4SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */ 16885add1ac4SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5), 16895add1ac4SLaurent Vivier status); /* A7 */ 16905add1ac4SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1), 16915add1ac4SLaurent Vivier status); /* A6 */ 16925add1ac4SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */ 16935add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */ 16945add1ac4SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 16955add1ac4SLaurent Vivier make_float64(0xBE5AE6452A118AE4), status), 16965add1ac4SLaurent Vivier status); /* A5+T*A7 */ 16975add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 16985add1ac4SLaurent Vivier make_float64(0x3EC71DE3A5341531), status), 16995add1ac4SLaurent Vivier status); /* A4+T*A6 */ 17005add1ac4SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */ 17015add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */ 17025add1ac4SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 17035add1ac4SLaurent Vivier make_float64(0xBF2A01A01A018B59), status), 17045add1ac4SLaurent Vivier status); /* A3+T(A5+TA7) */ 17055add1ac4SLaurent Vivier fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF)); 17065add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */ 17075add1ac4SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */ 17085add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */ 17095add1ac4SLaurent Vivier fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99)); 17105add1ac4SLaurent Vivier fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */ 17115add1ac4SLaurent Vivier fp1 = floatx80_add(fp1, fp2, 17125add1ac4SLaurent Vivier status); /* [A1+T(A3+T(A5+TA7))]+ 17135add1ac4SLaurent Vivier * [S(A2+T(A4+TA6))] 17145add1ac4SLaurent Vivier */ 17155add1ac4SLaurent Vivier 17165add1ac4SLaurent Vivier x = packFloatx80(xSign, xExp, xSig); 17175add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, x, status); /* R'*S */ 17185add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */ 17195add1ac4SLaurent Vivier 17205add1ac4SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 17215add1ac4SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 17225add1ac4SLaurent Vivier 17235add1ac4SLaurent Vivier a = floatx80_add(fp0, x, status); 17245add1ac4SLaurent Vivier 17255add1ac4SLaurent Vivier float_raise(float_flag_inexact, status); 17265add1ac4SLaurent Vivier 17275add1ac4SLaurent Vivier return a; 17285add1ac4SLaurent Vivier } 17295add1ac4SLaurent Vivier } 17305add1ac4SLaurent Vivier } 173168d0ed37SLaurent Vivier 1732*808d77bcSLucien Murray-Pitts /* 1733*808d77bcSLucien Murray-Pitts * Cosine 1734*808d77bcSLucien Murray-Pitts */ 173568d0ed37SLaurent Vivier 173668d0ed37SLaurent Vivier floatx80 floatx80_cos(floatx80 a, float_status *status) 173768d0ed37SLaurent Vivier { 173868d0ed37SLaurent Vivier flag aSign, xSign; 173968d0ed37SLaurent Vivier int32_t aExp, xExp; 174068d0ed37SLaurent Vivier uint64_t aSig, xSig; 174168d0ed37SLaurent Vivier 174268d0ed37SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 174368d0ed37SLaurent Vivier 174468d0ed37SLaurent Vivier int32_t compact, l, n, j; 174568d0ed37SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2; 174668d0ed37SLaurent Vivier float32 posneg1, twoto63; 17476361d298SLaurent Vivier flag endflag; 174868d0ed37SLaurent Vivier 174968d0ed37SLaurent Vivier aSig = extractFloatx80Frac(a); 175068d0ed37SLaurent Vivier aExp = extractFloatx80Exp(a); 175168d0ed37SLaurent Vivier aSign = extractFloatx80Sign(a); 175268d0ed37SLaurent Vivier 175368d0ed37SLaurent Vivier if (aExp == 0x7FFF) { 175468d0ed37SLaurent Vivier if ((uint64_t) (aSig << 1)) { 175568d0ed37SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 175668d0ed37SLaurent Vivier } 175768d0ed37SLaurent Vivier float_raise(float_flag_invalid, status); 175868d0ed37SLaurent Vivier return floatx80_default_nan(status); 175968d0ed37SLaurent Vivier } 176068d0ed37SLaurent Vivier 176168d0ed37SLaurent Vivier if (aExp == 0 && aSig == 0) { 176268d0ed37SLaurent Vivier return packFloatx80(0, one_exp, one_sig); 176368d0ed37SLaurent Vivier } 176468d0ed37SLaurent Vivier 176568d0ed37SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 176668d0ed37SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 176768d0ed37SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 176868d0ed37SLaurent Vivier status->floatx80_rounding_precision = 80; 176968d0ed37SLaurent Vivier 177068d0ed37SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 177168d0ed37SLaurent Vivier 177268d0ed37SLaurent Vivier fp0 = a; 177368d0ed37SLaurent Vivier 177468d0ed37SLaurent Vivier if (compact < 0x3FD78000 || compact > 0x4004BC7E) { 177568d0ed37SLaurent Vivier /* 2^(-40) > |X| > 15 PI */ 177668d0ed37SLaurent Vivier if (compact > 0x3FFF8000) { /* |X| >= 15 PI */ 177768d0ed37SLaurent Vivier /* REDUCEX */ 177868d0ed37SLaurent Vivier fp1 = packFloatx80(0, 0, 0); 177968d0ed37SLaurent Vivier if (compact == 0x7FFEFFFF) { 178068d0ed37SLaurent Vivier twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, 178168d0ed37SLaurent Vivier LIT64(0xC90FDAA200000000)); 178268d0ed37SLaurent Vivier twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, 178368d0ed37SLaurent Vivier LIT64(0x85A308D300000000)); 178468d0ed37SLaurent Vivier fp0 = floatx80_add(fp0, twopi1, status); 178568d0ed37SLaurent Vivier fp1 = fp0; 178668d0ed37SLaurent Vivier fp0 = floatx80_add(fp0, twopi2, status); 178768d0ed37SLaurent Vivier fp1 = floatx80_sub(fp1, fp0, status); 178868d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, twopi2, status); 178968d0ed37SLaurent Vivier } 179068d0ed37SLaurent Vivier loop: 179168d0ed37SLaurent Vivier xSign = extractFloatx80Sign(fp0); 179268d0ed37SLaurent Vivier xExp = extractFloatx80Exp(fp0); 179368d0ed37SLaurent Vivier xExp -= 0x3FFF; 179468d0ed37SLaurent Vivier if (xExp <= 28) { 179568d0ed37SLaurent Vivier l = 0; 179668d0ed37SLaurent Vivier endflag = 1; 179768d0ed37SLaurent Vivier } else { 179868d0ed37SLaurent Vivier l = xExp - 27; 179968d0ed37SLaurent Vivier endflag = 0; 180068d0ed37SLaurent Vivier } 180168d0ed37SLaurent Vivier invtwopi = packFloatx80(0, 0x3FFE - l, 180268d0ed37SLaurent Vivier LIT64(0xA2F9836E4E44152A)); /* INVTWOPI */ 180368d0ed37SLaurent Vivier twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); 180468d0ed37SLaurent Vivier twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); 180568d0ed37SLaurent Vivier 180668d0ed37SLaurent Vivier /* SIGN(INARG)*2^63 IN SGL */ 180768d0ed37SLaurent Vivier twoto63 = packFloat32(xSign, 0xBE, 0); 180868d0ed37SLaurent Vivier 180968d0ed37SLaurent Vivier fp2 = floatx80_mul(fp0, invtwopi, status); 181068d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), 181168d0ed37SLaurent Vivier status); /* THE FRACT PART OF FP2 IS ROUNDED */ 181268d0ed37SLaurent Vivier fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), 181368d0ed37SLaurent Vivier status); /* FP2 is N */ 181468d0ed37SLaurent Vivier fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */ 181568d0ed37SLaurent Vivier fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */ 181668d0ed37SLaurent Vivier fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */ 181768d0ed37SLaurent Vivier fp4 = floatx80_sub(fp4, fp3, status); /* W-P */ 181868d0ed37SLaurent Vivier fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */ 181968d0ed37SLaurent Vivier fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */ 182068d0ed37SLaurent Vivier fp3 = fp0; /* FP3 is A */ 182168d0ed37SLaurent Vivier fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */ 182268d0ed37SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */ 182368d0ed37SLaurent Vivier 182468d0ed37SLaurent Vivier if (endflag > 0) { 182568d0ed37SLaurent Vivier n = floatx80_to_int32(fp2, status); 182668d0ed37SLaurent Vivier goto sincont; 182768d0ed37SLaurent Vivier } 182868d0ed37SLaurent Vivier fp3 = floatx80_sub(fp3, fp0, status); /* A-R */ 182968d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */ 183068d0ed37SLaurent Vivier goto loop; 183168d0ed37SLaurent Vivier } else { 183268d0ed37SLaurent Vivier /* SINSM */ 183368d0ed37SLaurent Vivier fp0 = float32_to_floatx80(make_float32(0x3F800000), status); /* 1 */ 183468d0ed37SLaurent Vivier 183568d0ed37SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 183668d0ed37SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 183768d0ed37SLaurent Vivier 183868d0ed37SLaurent Vivier /* COSTINY */ 183968d0ed37SLaurent Vivier a = floatx80_sub(fp0, float32_to_floatx80( 184068d0ed37SLaurent Vivier make_float32(0x00800000), status), 184168d0ed37SLaurent Vivier status); 184268d0ed37SLaurent Vivier float_raise(float_flag_inexact, status); 184368d0ed37SLaurent Vivier 184468d0ed37SLaurent Vivier return a; 184568d0ed37SLaurent Vivier } 184668d0ed37SLaurent Vivier } else { 184768d0ed37SLaurent Vivier fp1 = floatx80_mul(fp0, float64_to_floatx80( 184868d0ed37SLaurent Vivier make_float64(0x3FE45F306DC9C883), status), 184968d0ed37SLaurent Vivier status); /* X*2/PI */ 185068d0ed37SLaurent Vivier 185168d0ed37SLaurent Vivier n = floatx80_to_int32(fp1, status); 185268d0ed37SLaurent Vivier j = 32 + n; 185368d0ed37SLaurent Vivier 185468d0ed37SLaurent Vivier fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */ 185568d0ed37SLaurent Vivier fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status), 185668d0ed37SLaurent Vivier status); /* FP0 IS R = (X-Y1)-Y2 */ 185768d0ed37SLaurent Vivier 185868d0ed37SLaurent Vivier sincont: 18596361d298SLaurent Vivier if ((n + 1) & 1) { 186068d0ed37SLaurent Vivier /* COSPOLY */ 186168d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */ 186268d0ed37SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */ 186368d0ed37SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3), 186468d0ed37SLaurent Vivier status); /* B8 */ 186568d0ed37SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19), 186668d0ed37SLaurent Vivier status); /* B7 */ 186768d0ed37SLaurent Vivier 186868d0ed37SLaurent Vivier xSign = extractFloatx80Sign(fp0); /* X IS S */ 186968d0ed37SLaurent Vivier xExp = extractFloatx80Exp(fp0); 187068d0ed37SLaurent Vivier xSig = extractFloatx80Frac(fp0); 187168d0ed37SLaurent Vivier 18726361d298SLaurent Vivier if (((n + 1) >> 1) & 1) { 187368d0ed37SLaurent Vivier xSign ^= 1; 187468d0ed37SLaurent Vivier posneg1 = make_float32(0xBF800000); /* -1 */ 187568d0ed37SLaurent Vivier } else { 187668d0ed37SLaurent Vivier xSign ^= 0; 187768d0ed37SLaurent Vivier posneg1 = make_float32(0x3F800000); /* 1 */ 187868d0ed37SLaurent Vivier } /* X IS NOW R'= SGN*R */ 187968d0ed37SLaurent Vivier 188068d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */ 188168d0ed37SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */ 188268d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 188368d0ed37SLaurent Vivier make_float64(0x3E21EED90612C972), status), 188468d0ed37SLaurent Vivier status); /* B6+TB8 */ 188568d0ed37SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 188668d0ed37SLaurent Vivier make_float64(0xBE927E4FB79D9FCF), status), 188768d0ed37SLaurent Vivier status); /* B5+TB7 */ 188868d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */ 188968d0ed37SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */ 189068d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 189168d0ed37SLaurent Vivier make_float64(0x3EFA01A01A01D423), status), 189268d0ed37SLaurent Vivier status); /* B4+T(B6+TB8) */ 189368d0ed37SLaurent Vivier fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438)); 189468d0ed37SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */ 189568d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */ 189668d0ed37SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */ 189768d0ed37SLaurent Vivier fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E)); 189868d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */ 189968d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, float32_to_floatx80( 190068d0ed37SLaurent Vivier make_float32(0xBF000000), status), 190168d0ed37SLaurent Vivier status); /* B1+T(B3+T(B5+TB7)) */ 190268d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */ 190368d0ed37SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); 190468d0ed37SLaurent Vivier /* [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))] */ 190568d0ed37SLaurent Vivier 190668d0ed37SLaurent Vivier x = packFloatx80(xSign, xExp, xSig); 190768d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, x, status); 190868d0ed37SLaurent Vivier 190968d0ed37SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 191068d0ed37SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 191168d0ed37SLaurent Vivier 191268d0ed37SLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status); 191368d0ed37SLaurent Vivier 191468d0ed37SLaurent Vivier float_raise(float_flag_inexact, status); 191568d0ed37SLaurent Vivier 191668d0ed37SLaurent Vivier return a; 191768d0ed37SLaurent Vivier } else { 191868d0ed37SLaurent Vivier /* SINPOLY */ 191968d0ed37SLaurent Vivier xSign = extractFloatx80Sign(fp0); /* X IS R */ 192068d0ed37SLaurent Vivier xExp = extractFloatx80Exp(fp0); 192168d0ed37SLaurent Vivier xSig = extractFloatx80Frac(fp0); 192268d0ed37SLaurent Vivier 19236361d298SLaurent Vivier xSign ^= ((n + 1) >> 1) & 1; /* X IS NOW R'= SGN*R */ 192468d0ed37SLaurent Vivier 192568d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */ 192668d0ed37SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */ 192768d0ed37SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5), 192868d0ed37SLaurent Vivier status); /* A7 */ 192968d0ed37SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1), 193068d0ed37SLaurent Vivier status); /* A6 */ 193168d0ed37SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */ 193268d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */ 193368d0ed37SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 193468d0ed37SLaurent Vivier make_float64(0xBE5AE6452A118AE4), status), 193568d0ed37SLaurent Vivier status); /* A5+T*A7 */ 193668d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 193768d0ed37SLaurent Vivier make_float64(0x3EC71DE3A5341531), status), 193868d0ed37SLaurent Vivier status); /* A4+T*A6 */ 193968d0ed37SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */ 194068d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */ 194168d0ed37SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 194268d0ed37SLaurent Vivier make_float64(0xBF2A01A01A018B59), status), 194368d0ed37SLaurent Vivier status); /* A3+T(A5+TA7) */ 194468d0ed37SLaurent Vivier fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF)); 194568d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */ 194668d0ed37SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */ 194768d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */ 194868d0ed37SLaurent Vivier fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99)); 194968d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */ 195068d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); 195168d0ed37SLaurent Vivier /* [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))] */ 195268d0ed37SLaurent Vivier 195368d0ed37SLaurent Vivier x = packFloatx80(xSign, xExp, xSig); 195468d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, x, status); /* R'*S */ 195568d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */ 195668d0ed37SLaurent Vivier 195768d0ed37SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 195868d0ed37SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 195968d0ed37SLaurent Vivier 196068d0ed37SLaurent Vivier a = floatx80_add(fp0, x, status); 196168d0ed37SLaurent Vivier 196268d0ed37SLaurent Vivier float_raise(float_flag_inexact, status); 196368d0ed37SLaurent Vivier 196468d0ed37SLaurent Vivier return a; 196568d0ed37SLaurent Vivier } 196668d0ed37SLaurent Vivier } 196768d0ed37SLaurent Vivier } 19688c992abcSLaurent Vivier 1969*808d77bcSLucien Murray-Pitts /* 1970*808d77bcSLucien Murray-Pitts * Arc tangent 1971*808d77bcSLucien Murray-Pitts */ 19728c992abcSLaurent Vivier 19738c992abcSLaurent Vivier floatx80 floatx80_atan(floatx80 a, float_status *status) 19748c992abcSLaurent Vivier { 19758c992abcSLaurent Vivier flag aSign; 19768c992abcSLaurent Vivier int32_t aExp; 19778c992abcSLaurent Vivier uint64_t aSig; 19788c992abcSLaurent Vivier 19798c992abcSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 19808c992abcSLaurent Vivier 19818c992abcSLaurent Vivier int32_t compact, tbl_index; 19828c992abcSLaurent Vivier floatx80 fp0, fp1, fp2, fp3, xsave; 19838c992abcSLaurent Vivier 19848c992abcSLaurent Vivier aSig = extractFloatx80Frac(a); 19858c992abcSLaurent Vivier aExp = extractFloatx80Exp(a); 19868c992abcSLaurent Vivier aSign = extractFloatx80Sign(a); 19878c992abcSLaurent Vivier 19888c992abcSLaurent Vivier if (aExp == 0x7FFF) { 19898c992abcSLaurent Vivier if ((uint64_t) (aSig << 1)) { 19908c992abcSLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 19918c992abcSLaurent Vivier } 19928c992abcSLaurent Vivier a = packFloatx80(aSign, piby2_exp, pi_sig); 19938c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 19948c992abcSLaurent Vivier return floatx80_move(a, status); 19958c992abcSLaurent Vivier } 19968c992abcSLaurent Vivier 19978c992abcSLaurent Vivier if (aExp == 0 && aSig == 0) { 19988c992abcSLaurent Vivier return packFloatx80(aSign, 0, 0); 19998c992abcSLaurent Vivier } 20008c992abcSLaurent Vivier 20018c992abcSLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 20028c992abcSLaurent Vivier 20038c992abcSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 20048c992abcSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 20058c992abcSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 20068c992abcSLaurent Vivier status->floatx80_rounding_precision = 80; 20078c992abcSLaurent Vivier 20088c992abcSLaurent Vivier if (compact < 0x3FFB8000 || compact > 0x4002FFFF) { 20098c992abcSLaurent Vivier /* |X| >= 16 or |X| < 1/16 */ 20108c992abcSLaurent Vivier if (compact > 0x3FFF8000) { /* |X| >= 16 */ 20118c992abcSLaurent Vivier if (compact > 0x40638000) { /* |X| > 2^(100) */ 20128c992abcSLaurent Vivier fp0 = packFloatx80(aSign, piby2_exp, pi_sig); 20138c992abcSLaurent Vivier fp1 = packFloatx80(aSign, 0x0001, one_sig); 20148c992abcSLaurent Vivier 20158c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 20168c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 20178c992abcSLaurent Vivier 20188c992abcSLaurent Vivier a = floatx80_sub(fp0, fp1, status); 20198c992abcSLaurent Vivier 20208c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 20218c992abcSLaurent Vivier 20228c992abcSLaurent Vivier return a; 20238c992abcSLaurent Vivier } else { 20248c992abcSLaurent Vivier fp0 = a; 20258c992abcSLaurent Vivier fp1 = packFloatx80(1, one_exp, one_sig); /* -1 */ 20268c992abcSLaurent Vivier fp1 = floatx80_div(fp1, fp0, status); /* X' = -1/X */ 20278c992abcSLaurent Vivier xsave = fp1; 20288c992abcSLaurent Vivier fp0 = floatx80_mul(fp1, fp1, status); /* Y = X'*X' */ 20298c992abcSLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */ 20308c992abcSLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBFB70BF398539E6A), 20318c992abcSLaurent Vivier status); /* C5 */ 20328c992abcSLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3FBC7187962D1D7D), 20338c992abcSLaurent Vivier status); /* C4 */ 20348c992abcSLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* Z*C5 */ 20358c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* Z*C4 */ 20368c992abcSLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 20378c992abcSLaurent Vivier make_float64(0xBFC24924827107B8), status), 20388c992abcSLaurent Vivier status); /* C3+Z*C5 */ 20398c992abcSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 20408c992abcSLaurent Vivier make_float64(0x3FC999999996263E), status), 20418c992abcSLaurent Vivier status); /* C2+Z*C4 */ 20428c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* Z*(C3+Z*C5) */ 20438c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* Y*(C2+Z*C4) */ 20448c992abcSLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 20458c992abcSLaurent Vivier make_float64(0xBFD5555555555536), status), 20468c992abcSLaurent Vivier status); /* C1+Z*(C3+Z*C5) */ 20478c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, xsave, status); /* X'*Y */ 20488c992abcSLaurent Vivier /* [Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)] */ 20498c992abcSLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); 20508c992abcSLaurent Vivier /* X'*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) ?? */ 20518c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); 20528c992abcSLaurent Vivier fp0 = floatx80_add(fp0, xsave, status); 20538c992abcSLaurent Vivier fp1 = packFloatx80(aSign, piby2_exp, pi_sig); 20548c992abcSLaurent Vivier 20558c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 20568c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 20578c992abcSLaurent Vivier 20588c992abcSLaurent Vivier a = floatx80_add(fp0, fp1, status); 20598c992abcSLaurent Vivier 20608c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 20618c992abcSLaurent Vivier 20628c992abcSLaurent Vivier return a; 20638c992abcSLaurent Vivier } 20648c992abcSLaurent Vivier } else { /* |X| < 1/16 */ 20658c992abcSLaurent Vivier if (compact < 0x3FD78000) { /* |X| < 2^(-40) */ 20668c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 20678c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 20688c992abcSLaurent Vivier 20698c992abcSLaurent Vivier a = floatx80_move(a, status); 20708c992abcSLaurent Vivier 20718c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 20728c992abcSLaurent Vivier 20738c992abcSLaurent Vivier return a; 20748c992abcSLaurent Vivier } else { 20758c992abcSLaurent Vivier fp0 = a; 20768c992abcSLaurent Vivier xsave = a; 20778c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* Y = X*X */ 20788c992abcSLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */ 20798c992abcSLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3FB344447F876989), 20808c992abcSLaurent Vivier status); /* B6 */ 20818c992abcSLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBFB744EE7FAF45DB), 20828c992abcSLaurent Vivier status); /* B5 */ 20838c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* Z*B6 */ 20848c992abcSLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* Z*B5 */ 20858c992abcSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 20868c992abcSLaurent Vivier make_float64(0x3FBC71C646940220), status), 20878c992abcSLaurent Vivier status); /* B4+Z*B6 */ 20888c992abcSLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 20898c992abcSLaurent Vivier make_float64(0xBFC24924921872F9), 20908c992abcSLaurent Vivier status), status); /* B3+Z*B5 */ 20918c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* Z*(B4+Z*B6) */ 20928c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* Z*(B3+Z*B5) */ 20938c992abcSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 20948c992abcSLaurent Vivier make_float64(0x3FC9999999998FA9), status), 20958c992abcSLaurent Vivier status); /* B2+Z*(B4+Z*B6) */ 20968c992abcSLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 20978c992abcSLaurent Vivier make_float64(0xBFD5555555555555), status), 20988c992abcSLaurent Vivier status); /* B1+Z*(B3+Z*B5) */ 20998c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* Y*(B2+Z*(B4+Z*B6)) */ 21008c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, xsave, status); /* X*Y */ 21018c992abcSLaurent Vivier /* [B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))] */ 21028c992abcSLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); 21038c992abcSLaurent Vivier /* X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) */ 21048c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); 21058c992abcSLaurent Vivier 21068c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 21078c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 21088c992abcSLaurent Vivier 21098c992abcSLaurent Vivier a = floatx80_add(fp0, xsave, status); 21108c992abcSLaurent Vivier 21118c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 21128c992abcSLaurent Vivier 21138c992abcSLaurent Vivier return a; 21148c992abcSLaurent Vivier } 21158c992abcSLaurent Vivier } 21168c992abcSLaurent Vivier } else { 21178c992abcSLaurent Vivier aSig &= LIT64(0xF800000000000000); 21188c992abcSLaurent Vivier aSig |= LIT64(0x0400000000000000); 21198c992abcSLaurent Vivier xsave = packFloatx80(aSign, aExp, aSig); /* F */ 21208c992abcSLaurent Vivier fp0 = a; 21218c992abcSLaurent Vivier fp1 = a; /* X */ 21228c992abcSLaurent Vivier fp2 = packFloatx80(0, one_exp, one_sig); /* 1 */ 21238c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, xsave, status); /* X*F */ 21248c992abcSLaurent Vivier fp0 = floatx80_sub(fp0, xsave, status); /* X-F */ 21258c992abcSLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); /* 1 + X*F */ 21268c992abcSLaurent Vivier fp0 = floatx80_div(fp0, fp1, status); /* U = (X-F)/(1+X*F) */ 21278c992abcSLaurent Vivier 21288c992abcSLaurent Vivier tbl_index = compact; 21298c992abcSLaurent Vivier 21308c992abcSLaurent Vivier tbl_index &= 0x7FFF0000; 21318c992abcSLaurent Vivier tbl_index -= 0x3FFB0000; 21328c992abcSLaurent Vivier tbl_index >>= 1; 21338c992abcSLaurent Vivier tbl_index += compact & 0x00007800; 21348c992abcSLaurent Vivier tbl_index >>= 11; 21358c992abcSLaurent Vivier 21368c992abcSLaurent Vivier fp3 = atan_tbl[tbl_index]; 21378c992abcSLaurent Vivier 21388c992abcSLaurent Vivier fp3.high |= aSign ? 0x8000 : 0; /* ATAN(F) */ 21398c992abcSLaurent Vivier 21408c992abcSLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* V = U*U */ 21418c992abcSLaurent Vivier fp2 = float64_to_floatx80(make_float64(0xBFF6687E314987D8), 21428c992abcSLaurent Vivier status); /* A3 */ 21438c992abcSLaurent Vivier fp2 = floatx80_add(fp2, fp1, status); /* A3+V */ 21448c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* V*(A3+V) */ 21458c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); /* U*V */ 21468c992abcSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 21478c992abcSLaurent Vivier make_float64(0x4002AC6934A26DB3), status), 21488c992abcSLaurent Vivier status); /* A2+V*(A3+V) */ 21498c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 21508c992abcSLaurent Vivier make_float64(0xBFC2476F4E1DA28E), status), 21518c992abcSLaurent Vivier status); /* A1+U*V */ 21528c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, fp2, status); /* A1*U*V*(A2+V*(A3+V)) */ 21538c992abcSLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* ATAN(U) */ 21548c992abcSLaurent Vivier 21558c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 21568c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 21578c992abcSLaurent Vivier 21588c992abcSLaurent Vivier a = floatx80_add(fp0, fp3, status); /* ATAN(X) */ 21598c992abcSLaurent Vivier 21608c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 21618c992abcSLaurent Vivier 21628c992abcSLaurent Vivier return a; 21638c992abcSLaurent Vivier } 21648c992abcSLaurent Vivier } 2165bc20b34eSLaurent Vivier 2166*808d77bcSLucien Murray-Pitts /* 2167*808d77bcSLucien Murray-Pitts * Arc sine 2168*808d77bcSLucien Murray-Pitts */ 2169bc20b34eSLaurent Vivier 2170bc20b34eSLaurent Vivier floatx80 floatx80_asin(floatx80 a, float_status *status) 2171bc20b34eSLaurent Vivier { 2172bc20b34eSLaurent Vivier flag aSign; 2173bc20b34eSLaurent Vivier int32_t aExp; 2174bc20b34eSLaurent Vivier uint64_t aSig; 2175bc20b34eSLaurent Vivier 2176bc20b34eSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 2177bc20b34eSLaurent Vivier 2178bc20b34eSLaurent Vivier int32_t compact; 2179bc20b34eSLaurent Vivier floatx80 fp0, fp1, fp2, one; 2180bc20b34eSLaurent Vivier 2181bc20b34eSLaurent Vivier aSig = extractFloatx80Frac(a); 2182bc20b34eSLaurent Vivier aExp = extractFloatx80Exp(a); 2183bc20b34eSLaurent Vivier aSign = extractFloatx80Sign(a); 2184bc20b34eSLaurent Vivier 2185bc20b34eSLaurent Vivier if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) { 2186bc20b34eSLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 2187bc20b34eSLaurent Vivier } 2188bc20b34eSLaurent Vivier 2189bc20b34eSLaurent Vivier if (aExp == 0 && aSig == 0) { 2190bc20b34eSLaurent Vivier return packFloatx80(aSign, 0, 0); 2191bc20b34eSLaurent Vivier } 2192bc20b34eSLaurent Vivier 2193bc20b34eSLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 2194bc20b34eSLaurent Vivier 2195bc20b34eSLaurent Vivier if (compact >= 0x3FFF8000) { /* |X| >= 1 */ 2196bc20b34eSLaurent Vivier if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */ 2197bc20b34eSLaurent Vivier float_raise(float_flag_inexact, status); 2198bc20b34eSLaurent Vivier a = packFloatx80(aSign, piby2_exp, pi_sig); 2199bc20b34eSLaurent Vivier return floatx80_move(a, status); 2200bc20b34eSLaurent Vivier } else { /* |X| > 1 */ 2201bc20b34eSLaurent Vivier float_raise(float_flag_invalid, status); 2202bc20b34eSLaurent Vivier return floatx80_default_nan(status); 2203bc20b34eSLaurent Vivier } 2204bc20b34eSLaurent Vivier 2205bc20b34eSLaurent Vivier } /* |X| < 1 */ 2206bc20b34eSLaurent Vivier 2207bc20b34eSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 2208bc20b34eSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 2209bc20b34eSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 2210bc20b34eSLaurent Vivier status->floatx80_rounding_precision = 80; 2211bc20b34eSLaurent Vivier 2212bc20b34eSLaurent Vivier one = packFloatx80(0, one_exp, one_sig); 2213bc20b34eSLaurent Vivier fp0 = a; 2214bc20b34eSLaurent Vivier 2215bc20b34eSLaurent Vivier fp1 = floatx80_sub(one, fp0, status); /* 1 - X */ 2216bc20b34eSLaurent Vivier fp2 = floatx80_add(one, fp0, status); /* 1 + X */ 2217bc20b34eSLaurent Vivier fp1 = floatx80_mul(fp2, fp1, status); /* (1+X)*(1-X) */ 2218bc20b34eSLaurent Vivier fp1 = floatx80_sqrt(fp1, status); /* SQRT((1+X)*(1-X)) */ 2219bc20b34eSLaurent Vivier fp0 = floatx80_div(fp0, fp1, status); /* X/SQRT((1+X)*(1-X)) */ 2220bc20b34eSLaurent Vivier 2221bc20b34eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 2222bc20b34eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 2223bc20b34eSLaurent Vivier 2224bc20b34eSLaurent Vivier a = floatx80_atan(fp0, status); /* ATAN(X/SQRT((1+X)*(1-X))) */ 2225bc20b34eSLaurent Vivier 2226bc20b34eSLaurent Vivier float_raise(float_flag_inexact, status); 2227bc20b34eSLaurent Vivier 2228bc20b34eSLaurent Vivier return a; 2229bc20b34eSLaurent Vivier } 2230c84813b8SLaurent Vivier 2231*808d77bcSLucien Murray-Pitts /* 2232*808d77bcSLucien Murray-Pitts * Arc cosine 2233*808d77bcSLucien Murray-Pitts */ 2234c84813b8SLaurent Vivier 2235c84813b8SLaurent Vivier floatx80 floatx80_acos(floatx80 a, float_status *status) 2236c84813b8SLaurent Vivier { 2237c84813b8SLaurent Vivier flag aSign; 2238c84813b8SLaurent Vivier int32_t aExp; 2239c84813b8SLaurent Vivier uint64_t aSig; 2240c84813b8SLaurent Vivier 2241c84813b8SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 2242c84813b8SLaurent Vivier 2243c84813b8SLaurent Vivier int32_t compact; 2244c84813b8SLaurent Vivier floatx80 fp0, fp1, one; 2245c84813b8SLaurent Vivier 2246c84813b8SLaurent Vivier aSig = extractFloatx80Frac(a); 2247c84813b8SLaurent Vivier aExp = extractFloatx80Exp(a); 2248c84813b8SLaurent Vivier aSign = extractFloatx80Sign(a); 2249c84813b8SLaurent Vivier 2250c84813b8SLaurent Vivier if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) { 2251c84813b8SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 2252c84813b8SLaurent Vivier } 2253c84813b8SLaurent Vivier if (aExp == 0 && aSig == 0) { 2254c84813b8SLaurent Vivier float_raise(float_flag_inexact, status); 2255c84813b8SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 0, 2256c84813b8SLaurent Vivier piby2_exp, pi_sig, 0, status); 2257c84813b8SLaurent Vivier } 2258c84813b8SLaurent Vivier 2259c84813b8SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 2260c84813b8SLaurent Vivier 2261c84813b8SLaurent Vivier if (compact >= 0x3FFF8000) { /* |X| >= 1 */ 2262c84813b8SLaurent Vivier if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */ 2263c84813b8SLaurent Vivier if (aSign) { /* X == -1 */ 2264c84813b8SLaurent Vivier a = packFloatx80(0, pi_exp, pi_sig); 2265c84813b8SLaurent Vivier float_raise(float_flag_inexact, status); 2266c84813b8SLaurent Vivier return floatx80_move(a, status); 2267c84813b8SLaurent Vivier } else { /* X == +1 */ 2268c84813b8SLaurent Vivier return packFloatx80(0, 0, 0); 2269c84813b8SLaurent Vivier } 2270c84813b8SLaurent Vivier } else { /* |X| > 1 */ 2271c84813b8SLaurent Vivier float_raise(float_flag_invalid, status); 2272c84813b8SLaurent Vivier return floatx80_default_nan(status); 2273c84813b8SLaurent Vivier } 2274c84813b8SLaurent Vivier } /* |X| < 1 */ 2275c84813b8SLaurent Vivier 2276c84813b8SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 2277c84813b8SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 2278c84813b8SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 2279c84813b8SLaurent Vivier status->floatx80_rounding_precision = 80; 2280c84813b8SLaurent Vivier 2281c84813b8SLaurent Vivier one = packFloatx80(0, one_exp, one_sig); 2282c84813b8SLaurent Vivier fp0 = a; 2283c84813b8SLaurent Vivier 2284c84813b8SLaurent Vivier fp1 = floatx80_add(one, fp0, status); /* 1 + X */ 2285c84813b8SLaurent Vivier fp0 = floatx80_sub(one, fp0, status); /* 1 - X */ 2286c84813b8SLaurent Vivier fp0 = floatx80_div(fp0, fp1, status); /* (1-X)/(1+X) */ 2287c84813b8SLaurent Vivier fp0 = floatx80_sqrt(fp0, status); /* SQRT((1-X)/(1+X)) */ 2288c84813b8SLaurent Vivier fp0 = floatx80_atan(fp0, status); /* ATAN(SQRT((1-X)/(1+X))) */ 2289c84813b8SLaurent Vivier 2290c84813b8SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 2291c84813b8SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 2292c84813b8SLaurent Vivier 2293c84813b8SLaurent Vivier a = floatx80_add(fp0, fp0, status); /* 2 * ATAN(SQRT((1-X)/(1+X))) */ 2294c84813b8SLaurent Vivier 2295c84813b8SLaurent Vivier float_raise(float_flag_inexact, status); 2296c84813b8SLaurent Vivier 2297c84813b8SLaurent Vivier return a; 2298c84813b8SLaurent Vivier } 2299e3655afaSLaurent Vivier 2300*808d77bcSLucien Murray-Pitts /* 2301*808d77bcSLucien Murray-Pitts * Hyperbolic arc tangent 2302*808d77bcSLucien Murray-Pitts */ 2303e3655afaSLaurent Vivier 2304e3655afaSLaurent Vivier floatx80 floatx80_atanh(floatx80 a, float_status *status) 2305e3655afaSLaurent Vivier { 2306e3655afaSLaurent Vivier flag aSign; 2307e3655afaSLaurent Vivier int32_t aExp; 2308e3655afaSLaurent Vivier uint64_t aSig; 2309e3655afaSLaurent Vivier 2310e3655afaSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 2311e3655afaSLaurent Vivier 2312e3655afaSLaurent Vivier int32_t compact; 2313e3655afaSLaurent Vivier floatx80 fp0, fp1, fp2, one; 2314e3655afaSLaurent Vivier 2315e3655afaSLaurent Vivier aSig = extractFloatx80Frac(a); 2316e3655afaSLaurent Vivier aExp = extractFloatx80Exp(a); 2317e3655afaSLaurent Vivier aSign = extractFloatx80Sign(a); 2318e3655afaSLaurent Vivier 2319e3655afaSLaurent Vivier if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) { 2320e3655afaSLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 2321e3655afaSLaurent Vivier } 2322e3655afaSLaurent Vivier 2323e3655afaSLaurent Vivier if (aExp == 0 && aSig == 0) { 2324e3655afaSLaurent Vivier return packFloatx80(aSign, 0, 0); 2325e3655afaSLaurent Vivier } 2326e3655afaSLaurent Vivier 2327e3655afaSLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 2328e3655afaSLaurent Vivier 2329e3655afaSLaurent Vivier if (compact >= 0x3FFF8000) { /* |X| >= 1 */ 2330e3655afaSLaurent Vivier if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */ 2331e3655afaSLaurent Vivier float_raise(float_flag_divbyzero, status); 2332e3655afaSLaurent Vivier return packFloatx80(aSign, floatx80_infinity.high, 2333e3655afaSLaurent Vivier floatx80_infinity.low); 2334e3655afaSLaurent Vivier } else { /* |X| > 1 */ 2335e3655afaSLaurent Vivier float_raise(float_flag_invalid, status); 2336e3655afaSLaurent Vivier return floatx80_default_nan(status); 2337e3655afaSLaurent Vivier } 2338e3655afaSLaurent Vivier } /* |X| < 1 */ 2339e3655afaSLaurent Vivier 2340e3655afaSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 2341e3655afaSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 2342e3655afaSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 2343e3655afaSLaurent Vivier status->floatx80_rounding_precision = 80; 2344e3655afaSLaurent Vivier 2345e3655afaSLaurent Vivier one = packFloatx80(0, one_exp, one_sig); 2346e3655afaSLaurent Vivier fp2 = packFloatx80(aSign, 0x3FFE, one_sig); /* SIGN(X) * (1/2) */ 2347e3655afaSLaurent Vivier fp0 = packFloatx80(0, aExp, aSig); /* Y = |X| */ 2348e3655afaSLaurent Vivier fp1 = packFloatx80(1, aExp, aSig); /* -Y */ 2349e3655afaSLaurent Vivier fp0 = floatx80_add(fp0, fp0, status); /* 2Y */ 2350e3655afaSLaurent Vivier fp1 = floatx80_add(fp1, one, status); /* 1-Y */ 2351e3655afaSLaurent Vivier fp0 = floatx80_div(fp0, fp1, status); /* Z = 2Y/(1-Y) */ 2352e3655afaSLaurent Vivier fp0 = floatx80_lognp1(fp0, status); /* LOG1P(Z) */ 2353e3655afaSLaurent Vivier 2354e3655afaSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 2355e3655afaSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 2356e3655afaSLaurent Vivier 2357e3655afaSLaurent Vivier a = floatx80_mul(fp0, fp2, 2358e3655afaSLaurent Vivier status); /* ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z) */ 2359e3655afaSLaurent Vivier 2360e3655afaSLaurent Vivier float_raise(float_flag_inexact, status); 2361e3655afaSLaurent Vivier 2362e3655afaSLaurent Vivier return a; 2363e3655afaSLaurent Vivier } 23649937b029SLaurent Vivier 2365*808d77bcSLucien Murray-Pitts /* 2366*808d77bcSLucien Murray-Pitts * e to x minus 1 2367*808d77bcSLucien Murray-Pitts */ 23689937b029SLaurent Vivier 23699937b029SLaurent Vivier floatx80 floatx80_etoxm1(floatx80 a, float_status *status) 23709937b029SLaurent Vivier { 23719937b029SLaurent Vivier flag aSign; 23729937b029SLaurent Vivier int32_t aExp; 23739937b029SLaurent Vivier uint64_t aSig; 23749937b029SLaurent Vivier 23759937b029SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 23769937b029SLaurent Vivier 23779937b029SLaurent Vivier int32_t compact, n, j, m, m1; 23789937b029SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, l2, sc, onebysc; 23799937b029SLaurent Vivier 23809937b029SLaurent Vivier aSig = extractFloatx80Frac(a); 23819937b029SLaurent Vivier aExp = extractFloatx80Exp(a); 23829937b029SLaurent Vivier aSign = extractFloatx80Sign(a); 23839937b029SLaurent Vivier 23849937b029SLaurent Vivier if (aExp == 0x7FFF) { 23859937b029SLaurent Vivier if ((uint64_t) (aSig << 1)) { 23869937b029SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 23879937b029SLaurent Vivier } 23889937b029SLaurent Vivier if (aSign) { 23899937b029SLaurent Vivier return packFloatx80(aSign, one_exp, one_sig); 23909937b029SLaurent Vivier } 23919937b029SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 23929937b029SLaurent Vivier floatx80_infinity.low); 23939937b029SLaurent Vivier } 23949937b029SLaurent Vivier 23959937b029SLaurent Vivier if (aExp == 0 && aSig == 0) { 23969937b029SLaurent Vivier return packFloatx80(aSign, 0, 0); 23979937b029SLaurent Vivier } 23989937b029SLaurent Vivier 23999937b029SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 24009937b029SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 24019937b029SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 24029937b029SLaurent Vivier status->floatx80_rounding_precision = 80; 24039937b029SLaurent Vivier 24049937b029SLaurent Vivier if (aExp >= 0x3FFD) { /* |X| >= 1/4 */ 24059937b029SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 24069937b029SLaurent Vivier 24079937b029SLaurent Vivier if (compact <= 0x4004C215) { /* |X| <= 70 log2 */ 24089937b029SLaurent Vivier fp0 = a; 24099937b029SLaurent Vivier fp1 = a; 24109937b029SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 24119937b029SLaurent Vivier make_float32(0x42B8AA3B), status), 24129937b029SLaurent Vivier status); /* 64/log2 * X */ 24139937b029SLaurent Vivier n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */ 24149937b029SLaurent Vivier fp0 = int32_to_floatx80(n, status); 24159937b029SLaurent Vivier 24169937b029SLaurent Vivier j = n & 0x3F; /* J = N mod 64 */ 24179937b029SLaurent Vivier m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */ 24189937b029SLaurent Vivier if (n < 0 && j) { 2419*808d77bcSLucien Murray-Pitts /* 2420*808d77bcSLucien Murray-Pitts * arithmetic right shift is division and 24219937b029SLaurent Vivier * round towards minus infinity 24229937b029SLaurent Vivier */ 24239937b029SLaurent Vivier m--; 24249937b029SLaurent Vivier } 24259937b029SLaurent Vivier m1 = -m; 24269937b029SLaurent Vivier /*m += 0x3FFF; // biased exponent of 2^(M) */ 24279937b029SLaurent Vivier /*m1 += 0x3FFF; // biased exponent of -2^(-M) */ 24289937b029SLaurent Vivier 24299937b029SLaurent Vivier fp2 = fp0; /* N */ 24309937b029SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 24319937b029SLaurent Vivier make_float32(0xBC317218), status), 24329937b029SLaurent Vivier status); /* N * L1, L1 = lead(-log2/64) */ 24339937b029SLaurent Vivier l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6)); 24349937b029SLaurent Vivier fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */ 24359937b029SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */ 24369937b029SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* R */ 24379937b029SLaurent Vivier 24389937b029SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 24399937b029SLaurent Vivier fp2 = float32_to_floatx80(make_float32(0x3950097B), 24409937b029SLaurent Vivier status); /* A6 */ 24419937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A6 */ 24429937b029SLaurent Vivier fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3AB60B6A), 24439937b029SLaurent Vivier status), fp1, status); /* fp3 is S*A5 */ 24449937b029SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 24459937b029SLaurent Vivier make_float64(0x3F81111111174385), status), 24469937b029SLaurent Vivier status); /* fp2 IS A4+S*A6 */ 24479937b029SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 24489937b029SLaurent Vivier make_float64(0x3FA5555555554F5A), status), 24499937b029SLaurent Vivier status); /* fp3 is A3+S*A5 */ 24509937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* fp2 IS S*(A4+S*A6) */ 24519937b029SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* fp3 IS S*(A3+S*A5) */ 24529937b029SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 24539937b029SLaurent Vivier make_float64(0x3FC5555555555555), status), 24549937b029SLaurent Vivier status); /* fp2 IS A2+S*(A4+S*A6) */ 24559937b029SLaurent Vivier fp3 = floatx80_add(fp3, float32_to_floatx80( 24569937b029SLaurent Vivier make_float32(0x3F000000), status), 24579937b029SLaurent Vivier status); /* fp3 IS A1+S*(A3+S*A5) */ 24589937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, 24599937b029SLaurent Vivier status); /* fp2 IS S*(A2+S*(A4+S*A6)) */ 24609937b029SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, 24619937b029SLaurent Vivier status); /* fp1 IS S*(A1+S*(A3+S*A5)) */ 24629937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, 24639937b029SLaurent Vivier status); /* fp2 IS R*S*(A2+S*(A4+S*A6)) */ 24649937b029SLaurent Vivier fp0 = floatx80_add(fp0, fp1, 24659937b029SLaurent Vivier status); /* fp0 IS R+S*(A1+S*(A3+S*A5)) */ 24669937b029SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */ 24679937b029SLaurent Vivier 24689937b029SLaurent Vivier fp0 = floatx80_mul(fp0, exp_tbl[j], 24699937b029SLaurent Vivier status); /* 2^(J/64)*(Exp(R)-1) */ 24709937b029SLaurent Vivier 24719937b029SLaurent Vivier if (m >= 64) { 24729937b029SLaurent Vivier fp1 = float32_to_floatx80(exp_tbl2[j], status); 24739937b029SLaurent Vivier onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */ 24749937b029SLaurent Vivier fp1 = floatx80_add(fp1, onebysc, status); 24759937b029SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); 24769937b029SLaurent Vivier fp0 = floatx80_add(fp0, exp_tbl[j], status); 24779937b029SLaurent Vivier } else if (m < -3) { 24789937b029SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], 24799937b029SLaurent Vivier status), status); 24809937b029SLaurent Vivier fp0 = floatx80_add(fp0, exp_tbl[j], status); 24819937b029SLaurent Vivier onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */ 24829937b029SLaurent Vivier fp0 = floatx80_add(fp0, onebysc, status); 24839937b029SLaurent Vivier } else { /* -3 <= m <= 63 */ 24849937b029SLaurent Vivier fp1 = exp_tbl[j]; 24859937b029SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], 24869937b029SLaurent Vivier status), status); 24879937b029SLaurent Vivier onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */ 24889937b029SLaurent Vivier fp1 = floatx80_add(fp1, onebysc, status); 24899937b029SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); 24909937b029SLaurent Vivier } 24919937b029SLaurent Vivier 24929937b029SLaurent Vivier sc = packFloatx80(0, m + 0x3FFF, one_sig); 24939937b029SLaurent Vivier 24949937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 24959937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 24969937b029SLaurent Vivier 24979937b029SLaurent Vivier a = floatx80_mul(fp0, sc, status); 24989937b029SLaurent Vivier 24999937b029SLaurent Vivier float_raise(float_flag_inexact, status); 25009937b029SLaurent Vivier 25019937b029SLaurent Vivier return a; 25029937b029SLaurent Vivier } else { /* |X| > 70 log2 */ 25039937b029SLaurent Vivier if (aSign) { 25049937b029SLaurent Vivier fp0 = float32_to_floatx80(make_float32(0xBF800000), 25059937b029SLaurent Vivier status); /* -1 */ 25069937b029SLaurent Vivier 25079937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 25089937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 25099937b029SLaurent Vivier 25109937b029SLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80( 25119937b029SLaurent Vivier make_float32(0x00800000), status), 25129937b029SLaurent Vivier status); /* -1 + 2^(-126) */ 25139937b029SLaurent Vivier 25149937b029SLaurent Vivier float_raise(float_flag_inexact, status); 25159937b029SLaurent Vivier 25169937b029SLaurent Vivier return a; 25179937b029SLaurent Vivier } else { 25189937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 25199937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 25209937b029SLaurent Vivier 25219937b029SLaurent Vivier return floatx80_etox(a, status); 25229937b029SLaurent Vivier } 25239937b029SLaurent Vivier } 25249937b029SLaurent Vivier } else { /* |X| < 1/4 */ 25259937b029SLaurent Vivier if (aExp >= 0x3FBE) { 25269937b029SLaurent Vivier fp0 = a; 25279937b029SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* S = X*X */ 25289937b029SLaurent Vivier fp1 = float32_to_floatx80(make_float32(0x2F30CAA8), 25299937b029SLaurent Vivier status); /* B12 */ 25309937b029SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); /* S * B12 */ 25319937b029SLaurent Vivier fp2 = float32_to_floatx80(make_float32(0x310F8290), 25329937b029SLaurent Vivier status); /* B11 */ 25339937b029SLaurent Vivier fp1 = floatx80_add(fp1, float32_to_floatx80( 25349937b029SLaurent Vivier make_float32(0x32D73220), status), 25359937b029SLaurent Vivier status); /* B10 */ 25369937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); 25379937b029SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); 25389937b029SLaurent Vivier fp2 = floatx80_add(fp2, float32_to_floatx80( 25399937b029SLaurent Vivier make_float32(0x3493F281), status), 25409937b029SLaurent Vivier status); /* B9 */ 25419937b029SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 25429937b029SLaurent Vivier make_float64(0x3EC71DE3A5774682), status), 25439937b029SLaurent Vivier status); /* B8 */ 25449937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); 25459937b029SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); 25469937b029SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 25479937b029SLaurent Vivier make_float64(0x3EFA01A019D7CB68), status), 25489937b029SLaurent Vivier status); /* B7 */ 25499937b029SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 25509937b029SLaurent Vivier make_float64(0x3F2A01A01A019DF3), status), 25519937b029SLaurent Vivier status); /* B6 */ 25529937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); 25539937b029SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); 25549937b029SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 25559937b029SLaurent Vivier make_float64(0x3F56C16C16C170E2), status), 25569937b029SLaurent Vivier status); /* B5 */ 25579937b029SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 25589937b029SLaurent Vivier make_float64(0x3F81111111111111), status), 25599937b029SLaurent Vivier status); /* B4 */ 25609937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); 25619937b029SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); 25629937b029SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 25639937b029SLaurent Vivier make_float64(0x3FA5555555555555), status), 25649937b029SLaurent Vivier status); /* B3 */ 25659937b029SLaurent Vivier fp3 = packFloatx80(0, 0x3FFC, LIT64(0xAAAAAAAAAAAAAAAB)); 25669937b029SLaurent Vivier fp1 = floatx80_add(fp1, fp3, status); /* B2 */ 25679937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); 25689937b029SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); 25699937b029SLaurent Vivier 25709937b029SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); 25719937b029SLaurent Vivier fp1 = floatx80_mul(fp1, a, status); 25729937b029SLaurent Vivier 25739937b029SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 25749937b029SLaurent Vivier make_float32(0x3F000000), status), 25759937b029SLaurent Vivier status); /* S*B1 */ 25769937b029SLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); /* Q */ 25779937b029SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* S*B1+Q */ 25789937b029SLaurent Vivier 25799937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 25809937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 25819937b029SLaurent Vivier 25829937b029SLaurent Vivier a = floatx80_add(fp0, a, status); 25839937b029SLaurent Vivier 25849937b029SLaurent Vivier float_raise(float_flag_inexact, status); 25859937b029SLaurent Vivier 25869937b029SLaurent Vivier return a; 25879937b029SLaurent Vivier } else { /* |X| < 2^(-65) */ 25889937b029SLaurent Vivier sc = packFloatx80(1, 1, one_sig); 25899937b029SLaurent Vivier fp0 = a; 25909937b029SLaurent Vivier 25919937b029SLaurent Vivier if (aExp < 0x0033) { /* |X| < 2^(-16382) */ 25929937b029SLaurent Vivier fp0 = floatx80_mul(fp0, float64_to_floatx80( 25939937b029SLaurent Vivier make_float64(0x48B0000000000000), status), 25949937b029SLaurent Vivier status); 25959937b029SLaurent Vivier fp0 = floatx80_add(fp0, sc, status); 25969937b029SLaurent Vivier 25979937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 25989937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 25999937b029SLaurent Vivier 26009937b029SLaurent Vivier a = floatx80_mul(fp0, float64_to_floatx80( 26019937b029SLaurent Vivier make_float64(0x3730000000000000), status), 26029937b029SLaurent Vivier status); 26039937b029SLaurent Vivier } else { 26049937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 26059937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 26069937b029SLaurent Vivier 26079937b029SLaurent Vivier a = floatx80_add(fp0, sc, status); 26089937b029SLaurent Vivier } 26099937b029SLaurent Vivier 26109937b029SLaurent Vivier float_raise(float_flag_inexact, status); 26119937b029SLaurent Vivier 26129937b029SLaurent Vivier return a; 26139937b029SLaurent Vivier } 26149937b029SLaurent Vivier } 26159937b029SLaurent Vivier } 26169937b029SLaurent Vivier 2617*808d77bcSLucien Murray-Pitts /* 2618*808d77bcSLucien Murray-Pitts * Hyperbolic tangent 2619*808d77bcSLucien Murray-Pitts */ 26209937b029SLaurent Vivier 26219937b029SLaurent Vivier floatx80 floatx80_tanh(floatx80 a, float_status *status) 26229937b029SLaurent Vivier { 26239937b029SLaurent Vivier flag aSign, vSign; 26249937b029SLaurent Vivier int32_t aExp, vExp; 26259937b029SLaurent Vivier uint64_t aSig, vSig; 26269937b029SLaurent Vivier 26279937b029SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 26289937b029SLaurent Vivier 26299937b029SLaurent Vivier int32_t compact; 26309937b029SLaurent Vivier floatx80 fp0, fp1; 26319937b029SLaurent Vivier uint32_t sign; 26329937b029SLaurent Vivier 26339937b029SLaurent Vivier aSig = extractFloatx80Frac(a); 26349937b029SLaurent Vivier aExp = extractFloatx80Exp(a); 26359937b029SLaurent Vivier aSign = extractFloatx80Sign(a); 26369937b029SLaurent Vivier 26379937b029SLaurent Vivier if (aExp == 0x7FFF) { 26389937b029SLaurent Vivier if ((uint64_t) (aSig << 1)) { 26399937b029SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 26409937b029SLaurent Vivier } 26419937b029SLaurent Vivier return packFloatx80(aSign, one_exp, one_sig); 26429937b029SLaurent Vivier } 26439937b029SLaurent Vivier 26449937b029SLaurent Vivier if (aExp == 0 && aSig == 0) { 26459937b029SLaurent Vivier return packFloatx80(aSign, 0, 0); 26469937b029SLaurent Vivier } 26479937b029SLaurent Vivier 26489937b029SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 26499937b029SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 26509937b029SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 26519937b029SLaurent Vivier status->floatx80_rounding_precision = 80; 26529937b029SLaurent Vivier 26539937b029SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 26549937b029SLaurent Vivier 26559937b029SLaurent Vivier if (compact < 0x3FD78000 || compact > 0x3FFFDDCE) { 26569937b029SLaurent Vivier /* TANHBORS */ 26579937b029SLaurent Vivier if (compact < 0x3FFF8000) { 26589937b029SLaurent Vivier /* TANHSM */ 26599937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 26609937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 26619937b029SLaurent Vivier 26629937b029SLaurent Vivier a = floatx80_move(a, status); 26639937b029SLaurent Vivier 26649937b029SLaurent Vivier float_raise(float_flag_inexact, status); 26659937b029SLaurent Vivier 26669937b029SLaurent Vivier return a; 26679937b029SLaurent Vivier } else { 26689937b029SLaurent Vivier if (compact > 0x40048AA1) { 26699937b029SLaurent Vivier /* TANHHUGE */ 26709937b029SLaurent Vivier sign = 0x3F800000; 26719937b029SLaurent Vivier sign |= aSign ? 0x80000000 : 0x00000000; 26729937b029SLaurent Vivier fp0 = float32_to_floatx80(make_float32(sign), status); 26739937b029SLaurent Vivier sign &= 0x80000000; 26749937b029SLaurent Vivier sign ^= 0x80800000; /* -SIGN(X)*EPS */ 26759937b029SLaurent Vivier 26769937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 26779937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 26789937b029SLaurent Vivier 26799937b029SLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80(make_float32(sign), 26809937b029SLaurent Vivier status), status); 26819937b029SLaurent Vivier 26829937b029SLaurent Vivier float_raise(float_flag_inexact, status); 26839937b029SLaurent Vivier 26849937b029SLaurent Vivier return a; 26859937b029SLaurent Vivier } else { 26869937b029SLaurent Vivier fp0 = packFloatx80(0, aExp + 1, aSig); /* Y = 2|X| */ 26879937b029SLaurent Vivier fp0 = floatx80_etox(fp0, status); /* FP0 IS EXP(Y) */ 26889937b029SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80( 26899937b029SLaurent Vivier make_float32(0x3F800000), 26909937b029SLaurent Vivier status), status); /* EXP(Y)+1 */ 26919937b029SLaurent Vivier sign = aSign ? 0x80000000 : 0x00000000; 26929937b029SLaurent Vivier fp1 = floatx80_div(float32_to_floatx80(make_float32( 26939937b029SLaurent Vivier sign ^ 0xC0000000), status), fp0, 26949937b029SLaurent Vivier status); /* -SIGN(X)*2 / [EXP(Y)+1] */ 26959937b029SLaurent Vivier fp0 = float32_to_floatx80(make_float32(sign | 0x3F800000), 26969937b029SLaurent Vivier status); /* SIGN */ 26979937b029SLaurent Vivier 26989937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 26999937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 27009937b029SLaurent Vivier 27019937b029SLaurent Vivier a = floatx80_add(fp1, fp0, status); 27029937b029SLaurent Vivier 27039937b029SLaurent Vivier float_raise(float_flag_inexact, status); 27049937b029SLaurent Vivier 27059937b029SLaurent Vivier return a; 27069937b029SLaurent Vivier } 27079937b029SLaurent Vivier } 27089937b029SLaurent Vivier } else { /* 2**(-40) < |X| < (5/2)LOG2 */ 27099937b029SLaurent Vivier fp0 = packFloatx80(0, aExp + 1, aSig); /* Y = 2|X| */ 27109937b029SLaurent Vivier fp0 = floatx80_etoxm1(fp0, status); /* FP0 IS Z = EXPM1(Y) */ 27119937b029SLaurent Vivier fp1 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x40000000), 27129937b029SLaurent Vivier status), 27139937b029SLaurent Vivier status); /* Z+2 */ 27149937b029SLaurent Vivier 27159937b029SLaurent Vivier vSign = extractFloatx80Sign(fp1); 27169937b029SLaurent Vivier vExp = extractFloatx80Exp(fp1); 27179937b029SLaurent Vivier vSig = extractFloatx80Frac(fp1); 27189937b029SLaurent Vivier 27199937b029SLaurent Vivier fp1 = packFloatx80(vSign ^ aSign, vExp, vSig); 27209937b029SLaurent Vivier 27219937b029SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 27229937b029SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 27239937b029SLaurent Vivier 27249937b029SLaurent Vivier a = floatx80_div(fp0, fp1, status); 27259937b029SLaurent Vivier 27269937b029SLaurent Vivier float_raise(float_flag_inexact, status); 27279937b029SLaurent Vivier 27289937b029SLaurent Vivier return a; 27299937b029SLaurent Vivier } 27309937b029SLaurent Vivier } 2731eee6b892SLaurent Vivier 2732*808d77bcSLucien Murray-Pitts /* 2733*808d77bcSLucien Murray-Pitts * Hyperbolic sine 2734*808d77bcSLucien Murray-Pitts */ 2735eee6b892SLaurent Vivier 2736eee6b892SLaurent Vivier floatx80 floatx80_sinh(floatx80 a, float_status *status) 2737eee6b892SLaurent Vivier { 2738eee6b892SLaurent Vivier flag aSign; 2739eee6b892SLaurent Vivier int32_t aExp; 2740eee6b892SLaurent Vivier uint64_t aSig; 2741eee6b892SLaurent Vivier 2742eee6b892SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 2743eee6b892SLaurent Vivier 2744eee6b892SLaurent Vivier int32_t compact; 2745eee6b892SLaurent Vivier floatx80 fp0, fp1, fp2; 2746eee6b892SLaurent Vivier float32 fact; 2747eee6b892SLaurent Vivier 2748eee6b892SLaurent Vivier aSig = extractFloatx80Frac(a); 2749eee6b892SLaurent Vivier aExp = extractFloatx80Exp(a); 2750eee6b892SLaurent Vivier aSign = extractFloatx80Sign(a); 2751eee6b892SLaurent Vivier 2752eee6b892SLaurent Vivier if (aExp == 0x7FFF) { 2753eee6b892SLaurent Vivier if ((uint64_t) (aSig << 1)) { 2754eee6b892SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 2755eee6b892SLaurent Vivier } 2756eee6b892SLaurent Vivier return packFloatx80(aSign, floatx80_infinity.high, 2757eee6b892SLaurent Vivier floatx80_infinity.low); 2758eee6b892SLaurent Vivier } 2759eee6b892SLaurent Vivier 2760eee6b892SLaurent Vivier if (aExp == 0 && aSig == 0) { 2761eee6b892SLaurent Vivier return packFloatx80(aSign, 0, 0); 2762eee6b892SLaurent Vivier } 2763eee6b892SLaurent Vivier 2764eee6b892SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 2765eee6b892SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 2766eee6b892SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 2767eee6b892SLaurent Vivier status->floatx80_rounding_precision = 80; 2768eee6b892SLaurent Vivier 2769eee6b892SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 2770eee6b892SLaurent Vivier 2771eee6b892SLaurent Vivier if (compact > 0x400CB167) { 2772eee6b892SLaurent Vivier /* SINHBIG */ 2773eee6b892SLaurent Vivier if (compact > 0x400CB2B3) { 2774eee6b892SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 2775eee6b892SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 2776eee6b892SLaurent Vivier 2777eee6b892SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 2778eee6b892SLaurent Vivier aSign, 0x8000, aSig, 0, status); 2779eee6b892SLaurent Vivier } else { 2780eee6b892SLaurent Vivier fp0 = floatx80_abs(a); /* Y = |X| */ 2781eee6b892SLaurent Vivier fp0 = floatx80_sub(fp0, float64_to_floatx80( 2782eee6b892SLaurent Vivier make_float64(0x40C62D38D3D64634), status), 2783eee6b892SLaurent Vivier status); /* (|X|-16381LOG2_LEAD) */ 2784eee6b892SLaurent Vivier fp0 = floatx80_sub(fp0, float64_to_floatx80( 2785eee6b892SLaurent Vivier make_float64(0x3D6F90AEB1E75CC7), status), 2786eee6b892SLaurent Vivier status); /* |X| - 16381 LOG2, ACCURATE */ 2787eee6b892SLaurent Vivier fp0 = floatx80_etox(fp0, status); 2788eee6b892SLaurent Vivier fp2 = packFloatx80(aSign, 0x7FFB, one_sig); 2789eee6b892SLaurent Vivier 2790eee6b892SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 2791eee6b892SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 2792eee6b892SLaurent Vivier 2793eee6b892SLaurent Vivier a = floatx80_mul(fp0, fp2, status); 2794eee6b892SLaurent Vivier 2795eee6b892SLaurent Vivier float_raise(float_flag_inexact, status); 2796eee6b892SLaurent Vivier 2797eee6b892SLaurent Vivier return a; 2798eee6b892SLaurent Vivier } 2799eee6b892SLaurent Vivier } else { /* |X| < 16380 LOG2 */ 2800eee6b892SLaurent Vivier fp0 = floatx80_abs(a); /* Y = |X| */ 2801eee6b892SLaurent Vivier fp0 = floatx80_etoxm1(fp0, status); /* FP0 IS Z = EXPM1(Y) */ 2802eee6b892SLaurent Vivier fp1 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000), 2803eee6b892SLaurent Vivier status), status); /* 1+Z */ 2804eee6b892SLaurent Vivier fp2 = fp0; 2805eee6b892SLaurent Vivier fp0 = floatx80_div(fp0, fp1, status); /* Z/(1+Z) */ 2806eee6b892SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); 2807eee6b892SLaurent Vivier 2808eee6b892SLaurent Vivier fact = packFloat32(aSign, 0x7E, 0); 2809eee6b892SLaurent Vivier 2810eee6b892SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 2811eee6b892SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 2812eee6b892SLaurent Vivier 2813eee6b892SLaurent Vivier a = floatx80_mul(fp0, float32_to_floatx80(fact, status), status); 2814eee6b892SLaurent Vivier 2815eee6b892SLaurent Vivier float_raise(float_flag_inexact, status); 2816eee6b892SLaurent Vivier 2817eee6b892SLaurent Vivier return a; 2818eee6b892SLaurent Vivier } 2819eee6b892SLaurent Vivier } 282002f9124eSLaurent Vivier 2821*808d77bcSLucien Murray-Pitts /* 2822*808d77bcSLucien Murray-Pitts * Hyperbolic cosine 2823*808d77bcSLucien Murray-Pitts */ 282402f9124eSLaurent Vivier 282502f9124eSLaurent Vivier floatx80 floatx80_cosh(floatx80 a, float_status *status) 282602f9124eSLaurent Vivier { 282702f9124eSLaurent Vivier int32_t aExp; 282802f9124eSLaurent Vivier uint64_t aSig; 282902f9124eSLaurent Vivier 283002f9124eSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 283102f9124eSLaurent Vivier 283202f9124eSLaurent Vivier int32_t compact; 283302f9124eSLaurent Vivier floatx80 fp0, fp1; 283402f9124eSLaurent Vivier 283502f9124eSLaurent Vivier aSig = extractFloatx80Frac(a); 283602f9124eSLaurent Vivier aExp = extractFloatx80Exp(a); 283702f9124eSLaurent Vivier 283802f9124eSLaurent Vivier if (aExp == 0x7FFF) { 283902f9124eSLaurent Vivier if ((uint64_t) (aSig << 1)) { 284002f9124eSLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 284102f9124eSLaurent Vivier } 284202f9124eSLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 284302f9124eSLaurent Vivier floatx80_infinity.low); 284402f9124eSLaurent Vivier } 284502f9124eSLaurent Vivier 284602f9124eSLaurent Vivier if (aExp == 0 && aSig == 0) { 284702f9124eSLaurent Vivier return packFloatx80(0, one_exp, one_sig); 284802f9124eSLaurent Vivier } 284902f9124eSLaurent Vivier 285002f9124eSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 285102f9124eSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 285202f9124eSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 285302f9124eSLaurent Vivier status->floatx80_rounding_precision = 80; 285402f9124eSLaurent Vivier 285502f9124eSLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 285602f9124eSLaurent Vivier 285702f9124eSLaurent Vivier if (compact > 0x400CB167) { 285802f9124eSLaurent Vivier if (compact > 0x400CB2B3) { 285902f9124eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 286002f9124eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 286102f9124eSLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 0, 286202f9124eSLaurent Vivier 0x8000, one_sig, 0, status); 286302f9124eSLaurent Vivier } else { 286402f9124eSLaurent Vivier fp0 = packFloatx80(0, aExp, aSig); 286502f9124eSLaurent Vivier fp0 = floatx80_sub(fp0, float64_to_floatx80( 286602f9124eSLaurent Vivier make_float64(0x40C62D38D3D64634), status), 286702f9124eSLaurent Vivier status); 286802f9124eSLaurent Vivier fp0 = floatx80_sub(fp0, float64_to_floatx80( 286902f9124eSLaurent Vivier make_float64(0x3D6F90AEB1E75CC7), status), 287002f9124eSLaurent Vivier status); 287102f9124eSLaurent Vivier fp0 = floatx80_etox(fp0, status); 287202f9124eSLaurent Vivier fp1 = packFloatx80(0, 0x7FFB, one_sig); 287302f9124eSLaurent Vivier 287402f9124eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 287502f9124eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 287602f9124eSLaurent Vivier 287702f9124eSLaurent Vivier a = floatx80_mul(fp0, fp1, status); 287802f9124eSLaurent Vivier 287902f9124eSLaurent Vivier float_raise(float_flag_inexact, status); 288002f9124eSLaurent Vivier 288102f9124eSLaurent Vivier return a; 288202f9124eSLaurent Vivier } 288302f9124eSLaurent Vivier } 288402f9124eSLaurent Vivier 288502f9124eSLaurent Vivier fp0 = packFloatx80(0, aExp, aSig); /* |X| */ 288602f9124eSLaurent Vivier fp0 = floatx80_etox(fp0, status); /* EXP(|X|) */ 288702f9124eSLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80(make_float32(0x3F000000), 288802f9124eSLaurent Vivier status), status); /* (1/2)*EXP(|X|) */ 288902f9124eSLaurent Vivier fp1 = float32_to_floatx80(make_float32(0x3E800000), status); /* 1/4 */ 289002f9124eSLaurent Vivier fp1 = floatx80_div(fp1, fp0, status); /* 1/(2*EXP(|X|)) */ 289102f9124eSLaurent Vivier 289202f9124eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 289302f9124eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 289402f9124eSLaurent Vivier 289502f9124eSLaurent Vivier a = floatx80_add(fp0, fp1, status); 289602f9124eSLaurent Vivier 289702f9124eSLaurent Vivier float_raise(float_flag_inexact, status); 289802f9124eSLaurent Vivier 289902f9124eSLaurent Vivier return a; 290002f9124eSLaurent Vivier } 2901