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