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