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 268c992abcSLaurent Vivier #define piby2_exp 0x3FFF 278c992abcSLaurent Vivier #define pi_sig LIT64(0xc90fdaa22168c235) 288c992abcSLaurent Vivier 290d379c17SLaurent Vivier static floatx80 propagateFloatx80NaNOneArg(floatx80 a, float_status *status) 300d379c17SLaurent Vivier { 310d379c17SLaurent Vivier if (floatx80_is_signaling_nan(a, status)) { 320d379c17SLaurent Vivier float_raise(float_flag_invalid, status); 330d379c17SLaurent Vivier } 340d379c17SLaurent Vivier 350d379c17SLaurent Vivier if (status->default_nan_mode) { 360d379c17SLaurent Vivier return floatx80_default_nan(status); 370d379c17SLaurent Vivier } 380d379c17SLaurent Vivier 390d379c17SLaurent Vivier return floatx80_maybe_silence_nan(a, status); 400d379c17SLaurent Vivier } 410d379c17SLaurent Vivier 42591596b7SLaurent Vivier /*---------------------------------------------------------------------------- 43591596b7SLaurent Vivier | Returns the modulo remainder of the extended double-precision floating-point 44591596b7SLaurent Vivier | value `a' with respect to the corresponding value `b'. 45591596b7SLaurent Vivier *----------------------------------------------------------------------------*/ 46591596b7SLaurent Vivier 47591596b7SLaurent Vivier floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status) 48591596b7SLaurent Vivier { 49591596b7SLaurent Vivier flag aSign, zSign; 50591596b7SLaurent Vivier int32_t aExp, bExp, expDiff; 51591596b7SLaurent Vivier uint64_t aSig0, aSig1, bSig; 52591596b7SLaurent Vivier uint64_t qTemp, term0, term1; 53591596b7SLaurent Vivier 54591596b7SLaurent Vivier aSig0 = extractFloatx80Frac(a); 55591596b7SLaurent Vivier aExp = extractFloatx80Exp(a); 56591596b7SLaurent Vivier aSign = extractFloatx80Sign(a); 57591596b7SLaurent Vivier bSig = extractFloatx80Frac(b); 58591596b7SLaurent Vivier bExp = extractFloatx80Exp(b); 59591596b7SLaurent Vivier 60591596b7SLaurent Vivier if (aExp == 0x7FFF) { 61591596b7SLaurent Vivier if ((uint64_t) (aSig0 << 1) 62591596b7SLaurent Vivier || ((bExp == 0x7FFF) && (uint64_t) (bSig << 1))) { 63591596b7SLaurent Vivier return propagateFloatx80NaN(a, b, status); 64591596b7SLaurent Vivier } 65591596b7SLaurent Vivier goto invalid; 66591596b7SLaurent Vivier } 67591596b7SLaurent Vivier if (bExp == 0x7FFF) { 68591596b7SLaurent Vivier if ((uint64_t) (bSig << 1)) { 69591596b7SLaurent Vivier return propagateFloatx80NaN(a, b, status); 70591596b7SLaurent Vivier } 71591596b7SLaurent Vivier return a; 72591596b7SLaurent Vivier } 73591596b7SLaurent Vivier if (bExp == 0) { 74591596b7SLaurent Vivier if (bSig == 0) { 75591596b7SLaurent Vivier invalid: 76591596b7SLaurent Vivier float_raise(float_flag_invalid, status); 77591596b7SLaurent Vivier return floatx80_default_nan(status); 78591596b7SLaurent Vivier } 79591596b7SLaurent Vivier normalizeFloatx80Subnormal(bSig, &bExp, &bSig); 80591596b7SLaurent Vivier } 81591596b7SLaurent Vivier if (aExp == 0) { 82591596b7SLaurent Vivier if ((uint64_t) (aSig0 << 1) == 0) { 83591596b7SLaurent Vivier return a; 84591596b7SLaurent Vivier } 85591596b7SLaurent Vivier normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0); 86591596b7SLaurent Vivier } 87591596b7SLaurent Vivier bSig |= LIT64(0x8000000000000000); 88591596b7SLaurent Vivier zSign = aSign; 89591596b7SLaurent Vivier expDiff = aExp - bExp; 90591596b7SLaurent Vivier aSig1 = 0; 91591596b7SLaurent Vivier if (expDiff < 0) { 92591596b7SLaurent Vivier return a; 93591596b7SLaurent Vivier } 94591596b7SLaurent Vivier qTemp = (bSig <= aSig0); 95591596b7SLaurent Vivier if (qTemp) { 96591596b7SLaurent Vivier aSig0 -= bSig; 97591596b7SLaurent Vivier } 98591596b7SLaurent Vivier expDiff -= 64; 99591596b7SLaurent Vivier while (0 < expDiff) { 100591596b7SLaurent Vivier qTemp = estimateDiv128To64(aSig0, aSig1, bSig); 101591596b7SLaurent Vivier qTemp = (2 < qTemp) ? qTemp - 2 : 0; 102591596b7SLaurent Vivier mul64To128(bSig, qTemp, &term0, &term1); 103591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 104591596b7SLaurent Vivier shortShift128Left(aSig0, aSig1, 62, &aSig0, &aSig1); 105591596b7SLaurent Vivier } 106591596b7SLaurent Vivier expDiff += 64; 107591596b7SLaurent Vivier if (0 < expDiff) { 108591596b7SLaurent Vivier qTemp = estimateDiv128To64(aSig0, aSig1, bSig); 109591596b7SLaurent Vivier qTemp = (2 < qTemp) ? qTemp - 2 : 0; 110591596b7SLaurent Vivier qTemp >>= 64 - expDiff; 111591596b7SLaurent Vivier mul64To128(bSig, qTemp << (64 - expDiff), &term0, &term1); 112591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 113591596b7SLaurent Vivier shortShift128Left(0, bSig, 64 - expDiff, &term0, &term1); 114591596b7SLaurent Vivier while (le128(term0, term1, aSig0, aSig1)) { 115591596b7SLaurent Vivier ++qTemp; 116591596b7SLaurent Vivier sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); 117591596b7SLaurent Vivier } 118591596b7SLaurent Vivier } 119591596b7SLaurent Vivier return 120591596b7SLaurent Vivier normalizeRoundAndPackFloatx80( 121591596b7SLaurent Vivier 80, zSign, bExp + expDiff, aSig0, aSig1, status); 122591596b7SLaurent Vivier } 1230d379c17SLaurent Vivier 1240d379c17SLaurent Vivier /*---------------------------------------------------------------------------- 1250d379c17SLaurent Vivier | Returns the mantissa of the extended double-precision floating-point 1260d379c17SLaurent Vivier | value `a'. 1270d379c17SLaurent Vivier *----------------------------------------------------------------------------*/ 1280d379c17SLaurent Vivier 1290d379c17SLaurent Vivier floatx80 floatx80_getman(floatx80 a, float_status *status) 1300d379c17SLaurent Vivier { 1310d379c17SLaurent Vivier flag aSign; 1320d379c17SLaurent Vivier int32_t aExp; 1330d379c17SLaurent Vivier uint64_t aSig; 1340d379c17SLaurent Vivier 1350d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 1360d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 1370d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 1380d379c17SLaurent Vivier 1390d379c17SLaurent Vivier if (aExp == 0x7FFF) { 1400d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 1410d379c17SLaurent Vivier return propagateFloatx80NaNOneArg(a , status); 1420d379c17SLaurent Vivier } 1430d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 1440d379c17SLaurent Vivier return floatx80_default_nan(status); 1450d379c17SLaurent Vivier } 1460d379c17SLaurent Vivier 1470d379c17SLaurent Vivier if (aExp == 0) { 1480d379c17SLaurent Vivier if (aSig == 0) { 1490d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 1500d379c17SLaurent Vivier } 1510d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 1520d379c17SLaurent Vivier } 1530d379c17SLaurent Vivier 1540d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign, 1550d379c17SLaurent Vivier 0x3FFF, aSig, 0, status); 1560d379c17SLaurent Vivier } 1570d379c17SLaurent Vivier 1580d379c17SLaurent Vivier /*---------------------------------------------------------------------------- 1590d379c17SLaurent Vivier | Returns the exponent of the extended double-precision floating-point 1600d379c17SLaurent Vivier | value `a' as an extended double-precision value. 1610d379c17SLaurent Vivier *----------------------------------------------------------------------------*/ 1620d379c17SLaurent Vivier 1630d379c17SLaurent Vivier floatx80 floatx80_getexp(floatx80 a, float_status *status) 1640d379c17SLaurent Vivier { 1650d379c17SLaurent Vivier flag aSign; 1660d379c17SLaurent Vivier int32_t aExp; 1670d379c17SLaurent Vivier uint64_t aSig; 1680d379c17SLaurent Vivier 1690d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 1700d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 1710d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 1720d379c17SLaurent Vivier 1730d379c17SLaurent Vivier if (aExp == 0x7FFF) { 1740d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 1750d379c17SLaurent Vivier return propagateFloatx80NaNOneArg(a , status); 1760d379c17SLaurent Vivier } 1770d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 1780d379c17SLaurent Vivier return floatx80_default_nan(status); 1790d379c17SLaurent Vivier } 1800d379c17SLaurent Vivier 1810d379c17SLaurent Vivier if (aExp == 0) { 1820d379c17SLaurent Vivier if (aSig == 0) { 1830d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 1840d379c17SLaurent Vivier } 1850d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 1860d379c17SLaurent Vivier } 1870d379c17SLaurent Vivier 1880d379c17SLaurent Vivier return int32_to_floatx80(aExp - 0x3FFF, status); 1890d379c17SLaurent Vivier } 1900d379c17SLaurent Vivier 1910d379c17SLaurent Vivier /*---------------------------------------------------------------------------- 1920d379c17SLaurent Vivier | Scales extended double-precision floating-point value in operand `a' by 1930d379c17SLaurent Vivier | value `b'. The function truncates the value in the second operand 'b' to 1940d379c17SLaurent Vivier | an integral value and adds that value to the exponent of the operand 'a'. 1950d379c17SLaurent Vivier | The operation performed according to the IEC/IEEE Standard for Binary 1960d379c17SLaurent Vivier | Floating-Point Arithmetic. 1970d379c17SLaurent Vivier *----------------------------------------------------------------------------*/ 1980d379c17SLaurent Vivier 1990d379c17SLaurent Vivier floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status) 2000d379c17SLaurent Vivier { 2010d379c17SLaurent Vivier flag aSign, bSign; 2020d379c17SLaurent Vivier int32_t aExp, bExp, shiftCount; 2030d379c17SLaurent Vivier uint64_t aSig, bSig; 2040d379c17SLaurent Vivier 2050d379c17SLaurent Vivier aSig = extractFloatx80Frac(a); 2060d379c17SLaurent Vivier aExp = extractFloatx80Exp(a); 2070d379c17SLaurent Vivier aSign = extractFloatx80Sign(a); 2080d379c17SLaurent Vivier bSig = extractFloatx80Frac(b); 2090d379c17SLaurent Vivier bExp = extractFloatx80Exp(b); 2100d379c17SLaurent Vivier bSign = extractFloatx80Sign(b); 2110d379c17SLaurent Vivier 2120d379c17SLaurent Vivier if (bExp == 0x7FFF) { 2130d379c17SLaurent Vivier if ((uint64_t) (bSig << 1) || 2140d379c17SLaurent Vivier ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) { 2150d379c17SLaurent Vivier return propagateFloatx80NaN(a, b, status); 2160d379c17SLaurent Vivier } 2170d379c17SLaurent Vivier float_raise(float_flag_invalid , status); 2180d379c17SLaurent Vivier return floatx80_default_nan(status); 2190d379c17SLaurent Vivier } 2200d379c17SLaurent Vivier if (aExp == 0x7FFF) { 2210d379c17SLaurent Vivier if ((uint64_t) (aSig << 1)) { 2220d379c17SLaurent Vivier return propagateFloatx80NaN(a, b, status); 2230d379c17SLaurent Vivier } 2240d379c17SLaurent Vivier return packFloatx80(aSign, floatx80_infinity.high, 2250d379c17SLaurent Vivier floatx80_infinity.low); 2260d379c17SLaurent Vivier } 2270d379c17SLaurent Vivier if (aExp == 0) { 2280d379c17SLaurent Vivier if (aSig == 0) { 2290d379c17SLaurent Vivier return packFloatx80(aSign, 0, 0); 2300d379c17SLaurent Vivier } 2310d379c17SLaurent Vivier if (bExp < 0x3FFF) { 2320d379c17SLaurent Vivier return a; 2330d379c17SLaurent Vivier } 2340d379c17SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 2350d379c17SLaurent Vivier } 2360d379c17SLaurent Vivier 2370d379c17SLaurent Vivier if (bExp < 0x3FFF) { 2380d379c17SLaurent Vivier return a; 2390d379c17SLaurent Vivier } 2400d379c17SLaurent Vivier 2410d379c17SLaurent Vivier if (0x400F < bExp) { 2420d379c17SLaurent Vivier aExp = bSign ? -0x6001 : 0xE000; 2430d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 2440d379c17SLaurent Vivier aSign, aExp, aSig, 0, status); 2450d379c17SLaurent Vivier } 2460d379c17SLaurent Vivier 2470d379c17SLaurent Vivier shiftCount = 0x403E - bExp; 2480d379c17SLaurent Vivier bSig >>= shiftCount; 2490d379c17SLaurent Vivier aExp = bSign ? (aExp - bSig) : (aExp + bSig); 2500d379c17SLaurent Vivier 2510d379c17SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 2520d379c17SLaurent Vivier aSign, aExp, aSig, 0, status); 2530d379c17SLaurent Vivier } 2549a069775SLaurent Vivier 2559a069775SLaurent Vivier floatx80 floatx80_move(floatx80 a, float_status *status) 2569a069775SLaurent Vivier { 2579a069775SLaurent Vivier flag aSign; 2589a069775SLaurent Vivier int32_t aExp; 2599a069775SLaurent Vivier uint64_t aSig; 2609a069775SLaurent Vivier 2619a069775SLaurent Vivier aSig = extractFloatx80Frac(a); 2629a069775SLaurent Vivier aExp = extractFloatx80Exp(a); 2639a069775SLaurent Vivier aSign = extractFloatx80Sign(a); 2649a069775SLaurent Vivier 2659a069775SLaurent Vivier if (aExp == 0x7FFF) { 2669a069775SLaurent Vivier if ((uint64_t)(aSig << 1)) { 2679a069775SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 2689a069775SLaurent Vivier } 2699a069775SLaurent Vivier return a; 2709a069775SLaurent Vivier } 2719a069775SLaurent Vivier if (aExp == 0) { 2729a069775SLaurent Vivier if (aSig == 0) { 2739a069775SLaurent Vivier return a; 2749a069775SLaurent Vivier } 2759a069775SLaurent Vivier normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision, 2769a069775SLaurent Vivier aSign, aExp, aSig, 0, status); 2779a069775SLaurent Vivier } 2789a069775SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign, 2799a069775SLaurent Vivier aExp, aSig, 0, status); 2809a069775SLaurent Vivier } 2814b5c65b8SLaurent Vivier 2824b5c65b8SLaurent Vivier /*---------------------------------------------------------------------------- 2834b5c65b8SLaurent Vivier | Algorithms for transcendental functions supported by MC68881 and MC68882 2844b5c65b8SLaurent Vivier | mathematical coprocessors. The functions are derived from FPSP library. 2854b5c65b8SLaurent Vivier *----------------------------------------------------------------------------*/ 2864b5c65b8SLaurent Vivier 2874b5c65b8SLaurent Vivier #define one_exp 0x3FFF 2884b5c65b8SLaurent Vivier #define one_sig LIT64(0x8000000000000000) 2894b5c65b8SLaurent Vivier 2904b5c65b8SLaurent Vivier /*---------------------------------------------------------------------------- 2914b5c65b8SLaurent Vivier | Function for compactifying extended double-precision floating point values. 2924b5c65b8SLaurent Vivier *----------------------------------------------------------------------------*/ 2934b5c65b8SLaurent Vivier 2944b5c65b8SLaurent Vivier static int32_t floatx80_make_compact(int32_t aExp, uint64_t aSig) 2954b5c65b8SLaurent Vivier { 2964b5c65b8SLaurent Vivier return (aExp << 16) | (aSig >> 48); 2974b5c65b8SLaurent Vivier } 2984b5c65b8SLaurent Vivier 2994b5c65b8SLaurent Vivier /*---------------------------------------------------------------------------- 3004b5c65b8SLaurent Vivier | Log base e of x plus 1 3014b5c65b8SLaurent Vivier *----------------------------------------------------------------------------*/ 3024b5c65b8SLaurent Vivier 3034b5c65b8SLaurent Vivier floatx80 floatx80_lognp1(floatx80 a, float_status *status) 3044b5c65b8SLaurent Vivier { 3054b5c65b8SLaurent Vivier flag aSign; 3064b5c65b8SLaurent Vivier int32_t aExp; 3074b5c65b8SLaurent Vivier uint64_t aSig, fSig; 3084b5c65b8SLaurent Vivier 3094b5c65b8SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 3104b5c65b8SLaurent Vivier 3114b5c65b8SLaurent Vivier int32_t compact, j, k; 3124b5c65b8SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu; 3134b5c65b8SLaurent Vivier 3144b5c65b8SLaurent Vivier aSig = extractFloatx80Frac(a); 3154b5c65b8SLaurent Vivier aExp = extractFloatx80Exp(a); 3164b5c65b8SLaurent Vivier aSign = extractFloatx80Sign(a); 3174b5c65b8SLaurent Vivier 3184b5c65b8SLaurent Vivier if (aExp == 0x7FFF) { 3194b5c65b8SLaurent Vivier if ((uint64_t) (aSig << 1)) { 3204b5c65b8SLaurent Vivier propagateFloatx80NaNOneArg(a, status); 3214b5c65b8SLaurent Vivier } 3224b5c65b8SLaurent Vivier if (aSign) { 3234b5c65b8SLaurent Vivier float_raise(float_flag_invalid, status); 3244b5c65b8SLaurent Vivier return floatx80_default_nan(status); 3254b5c65b8SLaurent Vivier } 3264b5c65b8SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, floatx80_infinity.low); 3274b5c65b8SLaurent Vivier } 3284b5c65b8SLaurent Vivier 3294b5c65b8SLaurent Vivier if (aExp == 0 && aSig == 0) { 3304b5c65b8SLaurent Vivier return packFloatx80(aSign, 0, 0); 3314b5c65b8SLaurent Vivier } 3324b5c65b8SLaurent Vivier 3334b5c65b8SLaurent Vivier if (aSign && aExp >= one_exp) { 3344b5c65b8SLaurent Vivier if (aExp == one_exp && aSig == one_sig) { 3354b5c65b8SLaurent Vivier float_raise(float_flag_divbyzero, status); 3364b5c65b8SLaurent Vivier packFloatx80(aSign, floatx80_infinity.high, floatx80_infinity.low); 3374b5c65b8SLaurent Vivier } 3384b5c65b8SLaurent Vivier float_raise(float_flag_invalid, status); 3394b5c65b8SLaurent Vivier return floatx80_default_nan(status); 3404b5c65b8SLaurent Vivier } 3414b5c65b8SLaurent Vivier 3424b5c65b8SLaurent Vivier if (aExp < 0x3f99 || (aExp == 0x3f99 && aSig == one_sig)) { 3434b5c65b8SLaurent Vivier /* <= min threshold */ 3444b5c65b8SLaurent Vivier float_raise(float_flag_inexact, status); 3454b5c65b8SLaurent Vivier return floatx80_move(a, status); 3464b5c65b8SLaurent Vivier } 3474b5c65b8SLaurent Vivier 3484b5c65b8SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 3494b5c65b8SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 3504b5c65b8SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 3514b5c65b8SLaurent Vivier status->floatx80_rounding_precision = 80; 3524b5c65b8SLaurent Vivier 3534b5c65b8SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 3544b5c65b8SLaurent Vivier 3554b5c65b8SLaurent Vivier fp0 = a; /* Z */ 3564b5c65b8SLaurent Vivier fp1 = a; 3574b5c65b8SLaurent Vivier 3584b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000), 3594b5c65b8SLaurent Vivier status), status); /* X = (1+Z) */ 3604b5c65b8SLaurent Vivier 3614b5c65b8SLaurent Vivier aExp = extractFloatx80Exp(fp0); 3624b5c65b8SLaurent Vivier aSig = extractFloatx80Frac(fp0); 3634b5c65b8SLaurent Vivier 3644b5c65b8SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 3654b5c65b8SLaurent Vivier 3664b5c65b8SLaurent Vivier if (compact < 0x3FFE8000 || compact > 0x3FFFC000) { 3674b5c65b8SLaurent Vivier /* |X| < 1/2 or |X| > 3/2 */ 3684b5c65b8SLaurent Vivier k = aExp - 0x3FFF; 3694b5c65b8SLaurent Vivier fp1 = int32_to_floatx80(k, status); 3704b5c65b8SLaurent Vivier 3714b5c65b8SLaurent Vivier fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); 3724b5c65b8SLaurent Vivier j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */ 3734b5c65b8SLaurent Vivier 3744b5c65b8SLaurent Vivier f = packFloatx80(0, 0x3FFF, fSig); /* F */ 3754b5c65b8SLaurent Vivier fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */ 3764b5c65b8SLaurent Vivier 3774b5c65b8SLaurent Vivier fp0 = floatx80_sub(fp0, f, status); /* Y-F */ 3784b5c65b8SLaurent Vivier 3794b5c65b8SLaurent Vivier lp1cont1: 3804b5c65b8SLaurent Vivier /* LP1CONT1 */ 3814b5c65b8SLaurent Vivier fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */ 3824b5c65b8SLaurent Vivier logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); 3834b5c65b8SLaurent Vivier klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */ 3844b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */ 3854b5c65b8SLaurent Vivier 3864b5c65b8SLaurent Vivier fp3 = fp2; 3874b5c65b8SLaurent Vivier fp1 = fp2; 3884b5c65b8SLaurent Vivier 3894b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 3904b5c65b8SLaurent Vivier make_float64(0x3FC2499AB5E4040B), status), 3914b5c65b8SLaurent Vivier status); /* V*A6 */ 3924b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, float64_to_floatx80( 3934b5c65b8SLaurent Vivier make_float64(0xBFC555B5848CB7DB), status), 3944b5c65b8SLaurent Vivier status); /* V*A5 */ 3954b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 3964b5c65b8SLaurent Vivier make_float64(0x3FC99999987D8730), status), 3974b5c65b8SLaurent Vivier status); /* A4+V*A6 */ 3984b5c65b8SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 3994b5c65b8SLaurent Vivier make_float64(0xBFCFFFFFFF6F7E97), status), 4004b5c65b8SLaurent Vivier status); /* A3+V*A5 */ 4014b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */ 4024b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */ 4034b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 4044b5c65b8SLaurent Vivier make_float64(0x3FD55555555555A4), status), 4054b5c65b8SLaurent Vivier status); /* A2+V*(A4+V*A6) */ 4064b5c65b8SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 4074b5c65b8SLaurent Vivier make_float64(0xBFE0000000000008), status), 4084b5c65b8SLaurent Vivier status); /* A1+V*(A3+V*A5) */ 4094b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */ 4104b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */ 4114b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */ 4124b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */ 4134b5c65b8SLaurent Vivier 4144b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, log_tbl[j + 1], 4154b5c65b8SLaurent Vivier status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */ 4164b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */ 4174b5c65b8SLaurent Vivier 4184b5c65b8SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 4194b5c65b8SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 4204b5c65b8SLaurent Vivier 4214b5c65b8SLaurent Vivier a = floatx80_add(fp0, klog2, status); 4224b5c65b8SLaurent Vivier 4234b5c65b8SLaurent Vivier float_raise(float_flag_inexact, status); 4244b5c65b8SLaurent Vivier 4254b5c65b8SLaurent Vivier return a; 4264b5c65b8SLaurent Vivier } else if (compact < 0x3FFEF07D || compact > 0x3FFF8841) { 4274b5c65b8SLaurent Vivier /* |X| < 1/16 or |X| > -1/16 */ 4284b5c65b8SLaurent Vivier /* LP1CARE */ 4294b5c65b8SLaurent Vivier fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); 4304b5c65b8SLaurent Vivier f = packFloatx80(0, 0x3FFF, fSig); /* F */ 4314b5c65b8SLaurent Vivier j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */ 4324b5c65b8SLaurent Vivier 4334b5c65b8SLaurent Vivier if (compact >= 0x3FFF8000) { /* 1+Z >= 1 */ 4344b5c65b8SLaurent Vivier /* KISZERO */ 4354b5c65b8SLaurent Vivier fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x3F800000), 4364b5c65b8SLaurent Vivier status), f, status); /* 1-F */ 4374b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (1-F)+Z */ 4384b5c65b8SLaurent Vivier fp1 = packFloatx80(0, 0, 0); /* K = 0 */ 4394b5c65b8SLaurent Vivier } else { 4404b5c65b8SLaurent Vivier /* KISNEG */ 4414b5c65b8SLaurent Vivier fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x40000000), 4424b5c65b8SLaurent Vivier status), f, status); /* 2-F */ 4434b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, fp1, status); /* 2Z */ 4444b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (2-F)+2Z */ 4454b5c65b8SLaurent Vivier fp1 = packFloatx80(1, one_exp, one_sig); /* K = -1 */ 4464b5c65b8SLaurent Vivier } 4474b5c65b8SLaurent Vivier goto lp1cont1; 4484b5c65b8SLaurent Vivier } else { 4494b5c65b8SLaurent Vivier /* LP1ONE16 */ 4504b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2Z */ 4514b5c65b8SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000), 4524b5c65b8SLaurent Vivier status), status); /* FP0 IS 1+X */ 4534b5c65b8SLaurent Vivier 4544b5c65b8SLaurent Vivier /* LP1CONT2 */ 4554b5c65b8SLaurent Vivier fp1 = floatx80_div(fp1, fp0, status); /* U */ 4564b5c65b8SLaurent Vivier saveu = fp1; 4574b5c65b8SLaurent Vivier fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */ 4584b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */ 4594b5c65b8SLaurent Vivier 4604b5c65b8SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6), 4614b5c65b8SLaurent Vivier status); /* B5 */ 4624b5c65b8SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0), 4634b5c65b8SLaurent Vivier status); /* B4 */ 4644b5c65b8SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */ 4654b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */ 4664b5c65b8SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 4674b5c65b8SLaurent Vivier make_float64(0x3F624924928BCCFF), status), 4684b5c65b8SLaurent Vivier status); /* B3+W*B5 */ 4694b5c65b8SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 4704b5c65b8SLaurent Vivier make_float64(0x3F899999999995EC), status), 4714b5c65b8SLaurent Vivier status); /* B2+W*B4 */ 4724b5c65b8SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */ 4734b5c65b8SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */ 4744b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 4754b5c65b8SLaurent Vivier make_float64(0x3FB5555555555555), status), 4764b5c65b8SLaurent Vivier status); /* B1+W*(B3+W*B5) */ 4774b5c65b8SLaurent Vivier 4784b5c65b8SLaurent Vivier fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */ 4794b5c65b8SLaurent Vivier fp1 = floatx80_add(fp1, fp2, 4804b5c65b8SLaurent Vivier status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */ 4814b5c65b8SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, 4824b5c65b8SLaurent Vivier status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */ 4834b5c65b8SLaurent Vivier 4844b5c65b8SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 4854b5c65b8SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 4864b5c65b8SLaurent Vivier 4874b5c65b8SLaurent Vivier a = floatx80_add(fp0, saveu, status); 4884b5c65b8SLaurent Vivier 4894b5c65b8SLaurent Vivier /*if (!floatx80_is_zero(a)) { */ 4904b5c65b8SLaurent Vivier float_raise(float_flag_inexact, status); 4914b5c65b8SLaurent Vivier /*} */ 4924b5c65b8SLaurent Vivier 4934b5c65b8SLaurent Vivier return a; 4944b5c65b8SLaurent Vivier } 4954b5c65b8SLaurent Vivier } 49650067bd1SLaurent Vivier 49750067bd1SLaurent Vivier /*---------------------------------------------------------------------------- 49850067bd1SLaurent Vivier | Log base e 49950067bd1SLaurent Vivier *----------------------------------------------------------------------------*/ 50050067bd1SLaurent Vivier 50150067bd1SLaurent Vivier floatx80 floatx80_logn(floatx80 a, float_status *status) 50250067bd1SLaurent Vivier { 50350067bd1SLaurent Vivier flag aSign; 50450067bd1SLaurent Vivier int32_t aExp; 50550067bd1SLaurent Vivier uint64_t aSig, fSig; 50650067bd1SLaurent Vivier 50750067bd1SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 50850067bd1SLaurent Vivier 50950067bd1SLaurent Vivier int32_t compact, j, k, adjk; 51050067bd1SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu; 51150067bd1SLaurent Vivier 51250067bd1SLaurent Vivier aSig = extractFloatx80Frac(a); 51350067bd1SLaurent Vivier aExp = extractFloatx80Exp(a); 51450067bd1SLaurent Vivier aSign = extractFloatx80Sign(a); 51550067bd1SLaurent Vivier 51650067bd1SLaurent Vivier if (aExp == 0x7FFF) { 51750067bd1SLaurent Vivier if ((uint64_t) (aSig << 1)) { 51850067bd1SLaurent Vivier propagateFloatx80NaNOneArg(a, status); 51950067bd1SLaurent Vivier } 52050067bd1SLaurent Vivier if (aSign == 0) { 52150067bd1SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 52250067bd1SLaurent Vivier floatx80_infinity.low); 52350067bd1SLaurent Vivier } 52450067bd1SLaurent Vivier } 52550067bd1SLaurent Vivier 52650067bd1SLaurent Vivier adjk = 0; 52750067bd1SLaurent Vivier 52850067bd1SLaurent Vivier if (aExp == 0) { 52950067bd1SLaurent Vivier if (aSig == 0) { /* zero */ 53050067bd1SLaurent Vivier float_raise(float_flag_divbyzero, status); 53150067bd1SLaurent Vivier return packFloatx80(1, floatx80_infinity.high, 53250067bd1SLaurent Vivier floatx80_infinity.low); 53350067bd1SLaurent Vivier } 53450067bd1SLaurent Vivier if ((aSig & one_sig) == 0) { /* denormal */ 53550067bd1SLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 53650067bd1SLaurent Vivier adjk = -100; 53750067bd1SLaurent Vivier aExp += 100; 53850067bd1SLaurent Vivier a = packFloatx80(aSign, aExp, aSig); 53950067bd1SLaurent Vivier } 54050067bd1SLaurent Vivier } 54150067bd1SLaurent Vivier 54250067bd1SLaurent Vivier if (aSign) { 54350067bd1SLaurent Vivier float_raise(float_flag_invalid, status); 54450067bd1SLaurent Vivier return floatx80_default_nan(status); 54550067bd1SLaurent Vivier } 54650067bd1SLaurent Vivier 54750067bd1SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 54850067bd1SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 54950067bd1SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 55050067bd1SLaurent Vivier status->floatx80_rounding_precision = 80; 55150067bd1SLaurent Vivier 55250067bd1SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 55350067bd1SLaurent Vivier 55450067bd1SLaurent Vivier if (compact < 0x3FFEF07D || compact > 0x3FFF8841) { 55550067bd1SLaurent Vivier /* |X| < 15/16 or |X| > 17/16 */ 55650067bd1SLaurent Vivier k = aExp - 0x3FFF; 55750067bd1SLaurent Vivier k += adjk; 55850067bd1SLaurent Vivier fp1 = int32_to_floatx80(k, status); 55950067bd1SLaurent Vivier 56050067bd1SLaurent Vivier fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); 56150067bd1SLaurent Vivier j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */ 56250067bd1SLaurent Vivier 56350067bd1SLaurent Vivier f = packFloatx80(0, 0x3FFF, fSig); /* F */ 56450067bd1SLaurent Vivier fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */ 56550067bd1SLaurent Vivier 56650067bd1SLaurent Vivier fp0 = floatx80_sub(fp0, f, status); /* Y-F */ 56750067bd1SLaurent Vivier 56850067bd1SLaurent Vivier /* LP1CONT1 */ 56950067bd1SLaurent Vivier fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */ 57050067bd1SLaurent Vivier logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); 57150067bd1SLaurent Vivier klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */ 57250067bd1SLaurent Vivier fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */ 57350067bd1SLaurent Vivier 57450067bd1SLaurent Vivier fp3 = fp2; 57550067bd1SLaurent Vivier fp1 = fp2; 57650067bd1SLaurent Vivier 57750067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 57850067bd1SLaurent Vivier make_float64(0x3FC2499AB5E4040B), status), 57950067bd1SLaurent Vivier status); /* V*A6 */ 58050067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, float64_to_floatx80( 58150067bd1SLaurent Vivier make_float64(0xBFC555B5848CB7DB), status), 58250067bd1SLaurent Vivier status); /* V*A5 */ 58350067bd1SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 58450067bd1SLaurent Vivier make_float64(0x3FC99999987D8730), status), 58550067bd1SLaurent Vivier status); /* A4+V*A6 */ 58650067bd1SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 58750067bd1SLaurent Vivier make_float64(0xBFCFFFFFFF6F7E97), status), 58850067bd1SLaurent Vivier status); /* A3+V*A5 */ 58950067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */ 59050067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */ 59150067bd1SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 59250067bd1SLaurent Vivier make_float64(0x3FD55555555555A4), status), 59350067bd1SLaurent Vivier status); /* A2+V*(A4+V*A6) */ 59450067bd1SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 59550067bd1SLaurent Vivier make_float64(0xBFE0000000000008), status), 59650067bd1SLaurent Vivier status); /* A1+V*(A3+V*A5) */ 59750067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */ 59850067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */ 59950067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */ 60050067bd1SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */ 60150067bd1SLaurent Vivier 60250067bd1SLaurent Vivier fp1 = floatx80_add(fp1, log_tbl[j + 1], 60350067bd1SLaurent Vivier status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */ 60450067bd1SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */ 60550067bd1SLaurent Vivier 60650067bd1SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 60750067bd1SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 60850067bd1SLaurent Vivier 60950067bd1SLaurent Vivier a = floatx80_add(fp0, klog2, status); 61050067bd1SLaurent Vivier 61150067bd1SLaurent Vivier float_raise(float_flag_inexact, status); 61250067bd1SLaurent Vivier 61350067bd1SLaurent Vivier return a; 61450067bd1SLaurent Vivier } else { /* |X-1| >= 1/16 */ 61550067bd1SLaurent Vivier fp0 = a; 61650067bd1SLaurent Vivier fp1 = a; 61750067bd1SLaurent Vivier fp1 = floatx80_sub(fp1, float32_to_floatx80(make_float32(0x3F800000), 61850067bd1SLaurent Vivier status), status); /* FP1 IS X-1 */ 61950067bd1SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000), 62050067bd1SLaurent Vivier status), status); /* FP0 IS X+1 */ 62150067bd1SLaurent Vivier fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2(X-1) */ 62250067bd1SLaurent Vivier 62350067bd1SLaurent Vivier /* LP1CONT2 */ 62450067bd1SLaurent Vivier fp1 = floatx80_div(fp1, fp0, status); /* U */ 62550067bd1SLaurent Vivier saveu = fp1; 62650067bd1SLaurent Vivier fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */ 62750067bd1SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */ 62850067bd1SLaurent Vivier 62950067bd1SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6), 63050067bd1SLaurent Vivier status); /* B5 */ 63150067bd1SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0), 63250067bd1SLaurent Vivier status); /* B4 */ 63350067bd1SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */ 63450067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */ 63550067bd1SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 63650067bd1SLaurent Vivier make_float64(0x3F624924928BCCFF), status), 63750067bd1SLaurent Vivier status); /* B3+W*B5 */ 63850067bd1SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 63950067bd1SLaurent Vivier make_float64(0x3F899999999995EC), status), 64050067bd1SLaurent Vivier status); /* B2+W*B4 */ 64150067bd1SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */ 64250067bd1SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */ 64350067bd1SLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 64450067bd1SLaurent Vivier make_float64(0x3FB5555555555555), status), 64550067bd1SLaurent Vivier status); /* B1+W*(B3+W*B5) */ 64650067bd1SLaurent Vivier 64750067bd1SLaurent Vivier fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */ 64850067bd1SLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */ 64950067bd1SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, 65050067bd1SLaurent Vivier status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */ 65150067bd1SLaurent Vivier 65250067bd1SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 65350067bd1SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 65450067bd1SLaurent Vivier 65550067bd1SLaurent Vivier a = floatx80_add(fp0, saveu, status); 65650067bd1SLaurent Vivier 65750067bd1SLaurent Vivier /*if (!floatx80_is_zero(a)) { */ 65850067bd1SLaurent Vivier float_raise(float_flag_inexact, status); 65950067bd1SLaurent Vivier /*} */ 66050067bd1SLaurent Vivier 66150067bd1SLaurent Vivier return a; 66250067bd1SLaurent Vivier } 66350067bd1SLaurent Vivier } 664248efb66SLaurent Vivier 665248efb66SLaurent Vivier /*---------------------------------------------------------------------------- 666248efb66SLaurent Vivier | Log base 10 667248efb66SLaurent Vivier *----------------------------------------------------------------------------*/ 668248efb66SLaurent Vivier 669248efb66SLaurent Vivier floatx80 floatx80_log10(floatx80 a, float_status *status) 670248efb66SLaurent Vivier { 671248efb66SLaurent Vivier flag aSign; 672248efb66SLaurent Vivier int32_t aExp; 673248efb66SLaurent Vivier uint64_t aSig; 674248efb66SLaurent Vivier 675248efb66SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 676248efb66SLaurent Vivier 677248efb66SLaurent Vivier floatx80 fp0, fp1; 678248efb66SLaurent Vivier 679248efb66SLaurent Vivier aSig = extractFloatx80Frac(a); 680248efb66SLaurent Vivier aExp = extractFloatx80Exp(a); 681248efb66SLaurent Vivier aSign = extractFloatx80Sign(a); 682248efb66SLaurent Vivier 683248efb66SLaurent Vivier if (aExp == 0x7FFF) { 684248efb66SLaurent Vivier if ((uint64_t) (aSig << 1)) { 685248efb66SLaurent Vivier propagateFloatx80NaNOneArg(a, status); 686248efb66SLaurent Vivier } 687248efb66SLaurent Vivier if (aSign == 0) { 688248efb66SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 689248efb66SLaurent Vivier floatx80_infinity.low); 690248efb66SLaurent Vivier } 691248efb66SLaurent Vivier } 692248efb66SLaurent Vivier 693248efb66SLaurent Vivier if (aExp == 0 && aSig == 0) { 694248efb66SLaurent Vivier float_raise(float_flag_divbyzero, status); 695248efb66SLaurent Vivier return packFloatx80(1, floatx80_infinity.high, 696248efb66SLaurent Vivier floatx80_infinity.low); 697248efb66SLaurent Vivier } 698248efb66SLaurent Vivier 699248efb66SLaurent Vivier if (aSign) { 700248efb66SLaurent Vivier float_raise(float_flag_invalid, status); 701248efb66SLaurent Vivier return floatx80_default_nan(status); 702248efb66SLaurent Vivier } 703248efb66SLaurent Vivier 704248efb66SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 705248efb66SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 706248efb66SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 707248efb66SLaurent Vivier status->floatx80_rounding_precision = 80; 708248efb66SLaurent Vivier 709248efb66SLaurent Vivier fp0 = floatx80_logn(a, status); 710248efb66SLaurent Vivier fp1 = packFloatx80(0, 0x3FFD, LIT64(0xDE5BD8A937287195)); /* INV_L10 */ 711248efb66SLaurent Vivier 712248efb66SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 713248efb66SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 714248efb66SLaurent Vivier 715248efb66SLaurent Vivier a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L10 */ 716248efb66SLaurent Vivier 717248efb66SLaurent Vivier float_raise(float_flag_inexact, status); 718248efb66SLaurent Vivier 719248efb66SLaurent Vivier return a; 720248efb66SLaurent Vivier } 72167b453edSLaurent Vivier 72267b453edSLaurent Vivier /*---------------------------------------------------------------------------- 72367b453edSLaurent Vivier | Log base 2 72467b453edSLaurent Vivier *----------------------------------------------------------------------------*/ 72567b453edSLaurent Vivier 72667b453edSLaurent Vivier floatx80 floatx80_log2(floatx80 a, float_status *status) 72767b453edSLaurent Vivier { 72867b453edSLaurent Vivier flag aSign; 72967b453edSLaurent Vivier int32_t aExp; 73067b453edSLaurent Vivier uint64_t aSig; 73167b453edSLaurent Vivier 73267b453edSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 73367b453edSLaurent Vivier 73467b453edSLaurent Vivier floatx80 fp0, fp1; 73567b453edSLaurent Vivier 73667b453edSLaurent Vivier aSig = extractFloatx80Frac(a); 73767b453edSLaurent Vivier aExp = extractFloatx80Exp(a); 73867b453edSLaurent Vivier aSign = extractFloatx80Sign(a); 73967b453edSLaurent Vivier 74067b453edSLaurent Vivier if (aExp == 0x7FFF) { 74167b453edSLaurent Vivier if ((uint64_t) (aSig << 1)) { 74267b453edSLaurent Vivier propagateFloatx80NaNOneArg(a, status); 74367b453edSLaurent Vivier } 74467b453edSLaurent Vivier if (aSign == 0) { 74567b453edSLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 74667b453edSLaurent Vivier floatx80_infinity.low); 74767b453edSLaurent Vivier } 74867b453edSLaurent Vivier } 74967b453edSLaurent Vivier 75067b453edSLaurent Vivier if (aExp == 0) { 75167b453edSLaurent Vivier if (aSig == 0) { 75267b453edSLaurent Vivier float_raise(float_flag_divbyzero, status); 75367b453edSLaurent Vivier return packFloatx80(1, floatx80_infinity.high, 75467b453edSLaurent Vivier floatx80_infinity.low); 75567b453edSLaurent Vivier } 75667b453edSLaurent Vivier normalizeFloatx80Subnormal(aSig, &aExp, &aSig); 75767b453edSLaurent Vivier } 75867b453edSLaurent Vivier 75967b453edSLaurent Vivier if (aSign) { 76067b453edSLaurent Vivier float_raise(float_flag_invalid, status); 76167b453edSLaurent Vivier return floatx80_default_nan(status); 76267b453edSLaurent Vivier } 76367b453edSLaurent Vivier 76467b453edSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 76567b453edSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 76667b453edSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 76767b453edSLaurent Vivier status->floatx80_rounding_precision = 80; 76867b453edSLaurent Vivier 76967b453edSLaurent Vivier if (aSig == one_sig) { /* X is 2^k */ 77067b453edSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 77167b453edSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 77267b453edSLaurent Vivier 77367b453edSLaurent Vivier a = int32_to_floatx80(aExp - 0x3FFF, status); 77467b453edSLaurent Vivier } else { 77567b453edSLaurent Vivier fp0 = floatx80_logn(a, status); 77667b453edSLaurent Vivier fp1 = packFloatx80(0, 0x3FFF, LIT64(0xB8AA3B295C17F0BC)); /* INV_L2 */ 77767b453edSLaurent Vivier 77867b453edSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 77967b453edSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 78067b453edSLaurent Vivier 78167b453edSLaurent Vivier a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L2 */ 78267b453edSLaurent Vivier } 78367b453edSLaurent Vivier 78467b453edSLaurent Vivier float_raise(float_flag_inexact, status); 78567b453edSLaurent Vivier 78667b453edSLaurent Vivier return a; 78767b453edSLaurent Vivier } 78840ad0873SLaurent Vivier 78940ad0873SLaurent Vivier /*---------------------------------------------------------------------------- 79040ad0873SLaurent Vivier | e to x 79140ad0873SLaurent Vivier *----------------------------------------------------------------------------*/ 79240ad0873SLaurent Vivier 79340ad0873SLaurent Vivier floatx80 floatx80_etox(floatx80 a, float_status *status) 79440ad0873SLaurent Vivier { 79540ad0873SLaurent Vivier flag aSign; 79640ad0873SLaurent Vivier int32_t aExp; 79740ad0873SLaurent Vivier uint64_t aSig; 79840ad0873SLaurent Vivier 79940ad0873SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 80040ad0873SLaurent Vivier 80140ad0873SLaurent Vivier int32_t compact, n, j, k, m, m1; 80240ad0873SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, l2, scale, adjscale; 80340ad0873SLaurent Vivier flag adjflag; 80440ad0873SLaurent Vivier 80540ad0873SLaurent Vivier aSig = extractFloatx80Frac(a); 80640ad0873SLaurent Vivier aExp = extractFloatx80Exp(a); 80740ad0873SLaurent Vivier aSign = extractFloatx80Sign(a); 80840ad0873SLaurent Vivier 80940ad0873SLaurent Vivier if (aExp == 0x7FFF) { 81040ad0873SLaurent Vivier if ((uint64_t) (aSig << 1)) { 81140ad0873SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 81240ad0873SLaurent Vivier } 81340ad0873SLaurent Vivier if (aSign) { 81440ad0873SLaurent Vivier return packFloatx80(0, 0, 0); 81540ad0873SLaurent Vivier } 81640ad0873SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 81740ad0873SLaurent Vivier floatx80_infinity.low); 81840ad0873SLaurent Vivier } 81940ad0873SLaurent Vivier 82040ad0873SLaurent Vivier if (aExp == 0 && aSig == 0) { 82140ad0873SLaurent Vivier return packFloatx80(0, one_exp, one_sig); 82240ad0873SLaurent Vivier } 82340ad0873SLaurent Vivier 82440ad0873SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 82540ad0873SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 82640ad0873SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 82740ad0873SLaurent Vivier status->floatx80_rounding_precision = 80; 82840ad0873SLaurent Vivier 82940ad0873SLaurent Vivier adjflag = 0; 83040ad0873SLaurent Vivier 83140ad0873SLaurent Vivier if (aExp >= 0x3FBE) { /* |X| >= 2^(-65) */ 83240ad0873SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 83340ad0873SLaurent Vivier 83440ad0873SLaurent Vivier if (compact < 0x400CB167) { /* |X| < 16380 log2 */ 83540ad0873SLaurent Vivier fp0 = a; 83640ad0873SLaurent Vivier fp1 = a; 83740ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 83840ad0873SLaurent Vivier make_float32(0x42B8AA3B), status), 83940ad0873SLaurent Vivier status); /* 64/log2 * X */ 84040ad0873SLaurent Vivier adjflag = 0; 84140ad0873SLaurent Vivier n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */ 84240ad0873SLaurent Vivier fp0 = int32_to_floatx80(n, status); 84340ad0873SLaurent Vivier 84440ad0873SLaurent Vivier j = n & 0x3F; /* J = N mod 64 */ 84540ad0873SLaurent Vivier m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */ 84640ad0873SLaurent Vivier if (n < 0 && j) { 84740ad0873SLaurent Vivier /* arithmetic right shift is division and 84840ad0873SLaurent Vivier * round towards minus infinity 84940ad0873SLaurent Vivier */ 85040ad0873SLaurent Vivier m--; 85140ad0873SLaurent Vivier } 85240ad0873SLaurent Vivier m += 0x3FFF; /* biased exponent of 2^(M) */ 85340ad0873SLaurent Vivier 85440ad0873SLaurent Vivier expcont1: 85540ad0873SLaurent Vivier fp2 = fp0; /* N */ 85640ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 85740ad0873SLaurent Vivier make_float32(0xBC317218), status), 85840ad0873SLaurent Vivier status); /* N * L1, L1 = lead(-log2/64) */ 85940ad0873SLaurent Vivier l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6)); 86040ad0873SLaurent Vivier fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */ 86140ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */ 86240ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* R */ 86340ad0873SLaurent Vivier 86440ad0873SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 86540ad0873SLaurent Vivier fp2 = float32_to_floatx80(make_float32(0x3AB60B70), 86640ad0873SLaurent Vivier status); /* A5 */ 86740ad0873SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A5 */ 86840ad0873SLaurent Vivier fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3C088895), 86940ad0873SLaurent Vivier status), fp1, 87040ad0873SLaurent Vivier status); /* fp3 is S*A4 */ 87140ad0873SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80(make_float64( 87240ad0873SLaurent Vivier 0x3FA5555555554431), status), 87340ad0873SLaurent Vivier status); /* fp2 is A3+S*A5 */ 87440ad0873SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80(make_float64( 87540ad0873SLaurent Vivier 0x3FC5555555554018), status), 87640ad0873SLaurent Vivier status); /* fp3 is A2+S*A4 */ 87740ad0873SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*(A3+S*A5) */ 87840ad0873SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* fp3 is S*(A2+S*A4) */ 87940ad0873SLaurent Vivier fp2 = floatx80_add(fp2, float32_to_floatx80( 88040ad0873SLaurent Vivier make_float32(0x3F000000), status), 88140ad0873SLaurent Vivier status); /* fp2 is A1+S*(A3+S*A5) */ 88240ad0873SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* fp3 IS R*S*(A2+S*A4) */ 88340ad0873SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, 88440ad0873SLaurent Vivier status); /* fp2 IS S*(A1+S*(A3+S*A5)) */ 88540ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp3, status); /* fp0 IS R+R*S*(A2+S*A4) */ 88640ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */ 88740ad0873SLaurent Vivier 88840ad0873SLaurent Vivier fp1 = exp_tbl[j]; 88940ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); /* 2^(J/64)*(Exp(R)-1) */ 89040ad0873SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status), 89140ad0873SLaurent Vivier status); /* accurate 2^(J/64) */ 89240ad0873SLaurent Vivier fp0 = floatx80_add(fp0, fp1, 89340ad0873SLaurent Vivier status); /* 2^(J/64) + 2^(J/64)*(Exp(R)-1) */ 89440ad0873SLaurent Vivier 89540ad0873SLaurent Vivier scale = packFloatx80(0, m, one_sig); 89640ad0873SLaurent Vivier if (adjflag) { 89740ad0873SLaurent Vivier adjscale = packFloatx80(0, m1, one_sig); 89840ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, adjscale, status); 89940ad0873SLaurent Vivier } 90040ad0873SLaurent Vivier 90140ad0873SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 90240ad0873SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 90340ad0873SLaurent Vivier 90440ad0873SLaurent Vivier a = floatx80_mul(fp0, scale, status); 90540ad0873SLaurent Vivier 90640ad0873SLaurent Vivier float_raise(float_flag_inexact, status); 90740ad0873SLaurent Vivier 90840ad0873SLaurent Vivier return a; 90940ad0873SLaurent Vivier } else { /* |X| >= 16380 log2 */ 91040ad0873SLaurent Vivier if (compact > 0x400CB27C) { /* |X| >= 16480 log2 */ 91140ad0873SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 91240ad0873SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 91340ad0873SLaurent Vivier if (aSign) { 91440ad0873SLaurent Vivier a = roundAndPackFloatx80( 91540ad0873SLaurent Vivier status->floatx80_rounding_precision, 91640ad0873SLaurent Vivier 0, -0x1000, aSig, 0, status); 91740ad0873SLaurent Vivier } else { 91840ad0873SLaurent Vivier a = roundAndPackFloatx80( 91940ad0873SLaurent Vivier status->floatx80_rounding_precision, 92040ad0873SLaurent Vivier 0, 0x8000, aSig, 0, status); 92140ad0873SLaurent Vivier } 92240ad0873SLaurent Vivier float_raise(float_flag_inexact, status); 92340ad0873SLaurent Vivier 92440ad0873SLaurent Vivier return a; 92540ad0873SLaurent Vivier } else { 92640ad0873SLaurent Vivier fp0 = a; 92740ad0873SLaurent Vivier fp1 = a; 92840ad0873SLaurent Vivier fp0 = floatx80_mul(fp0, float32_to_floatx80( 92940ad0873SLaurent Vivier make_float32(0x42B8AA3B), status), 93040ad0873SLaurent Vivier status); /* 64/log2 * X */ 93140ad0873SLaurent Vivier adjflag = 1; 93240ad0873SLaurent Vivier n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */ 93340ad0873SLaurent Vivier fp0 = int32_to_floatx80(n, status); 93440ad0873SLaurent Vivier 93540ad0873SLaurent Vivier j = n & 0x3F; /* J = N mod 64 */ 93640ad0873SLaurent Vivier /* NOTE: this is really arithmetic right shift by 6 */ 93740ad0873SLaurent Vivier k = n / 64; 93840ad0873SLaurent Vivier if (n < 0 && j) { 93940ad0873SLaurent Vivier /* arithmetic right shift is division and 94040ad0873SLaurent Vivier * round towards minus infinity 94140ad0873SLaurent Vivier */ 94240ad0873SLaurent Vivier k--; 94340ad0873SLaurent Vivier } 94440ad0873SLaurent Vivier /* NOTE: this is really arithmetic right shift by 1 */ 94540ad0873SLaurent Vivier m1 = k / 2; 94640ad0873SLaurent Vivier if (k < 0 && (k & 1)) { 94740ad0873SLaurent Vivier /* arithmetic right shift is division and 94840ad0873SLaurent Vivier * round towards minus infinity 94940ad0873SLaurent Vivier */ 95040ad0873SLaurent Vivier m1--; 95140ad0873SLaurent Vivier } 95240ad0873SLaurent Vivier m = k - m1; 95340ad0873SLaurent Vivier m1 += 0x3FFF; /* biased exponent of 2^(M1) */ 95440ad0873SLaurent Vivier m += 0x3FFF; /* biased exponent of 2^(M) */ 95540ad0873SLaurent Vivier 95640ad0873SLaurent Vivier goto expcont1; 95740ad0873SLaurent Vivier } 95840ad0873SLaurent Vivier } 95940ad0873SLaurent Vivier } else { /* |X| < 2^(-65) */ 96040ad0873SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 96140ad0873SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 96240ad0873SLaurent Vivier 96340ad0873SLaurent Vivier a = floatx80_add(a, float32_to_floatx80(make_float32(0x3F800000), 96440ad0873SLaurent Vivier status), status); /* 1 + X */ 96540ad0873SLaurent Vivier 96640ad0873SLaurent Vivier float_raise(float_flag_inexact, status); 96740ad0873SLaurent Vivier 96840ad0873SLaurent Vivier return a; 96940ad0873SLaurent Vivier } 97040ad0873SLaurent Vivier } 971068f1615SLaurent Vivier 972068f1615SLaurent Vivier /*---------------------------------------------------------------------------- 973068f1615SLaurent Vivier | 2 to x 974068f1615SLaurent Vivier *----------------------------------------------------------------------------*/ 975068f1615SLaurent Vivier 976068f1615SLaurent Vivier floatx80 floatx80_twotox(floatx80 a, float_status *status) 977068f1615SLaurent Vivier { 978068f1615SLaurent Vivier flag aSign; 979068f1615SLaurent Vivier int32_t aExp; 980068f1615SLaurent Vivier uint64_t aSig; 981068f1615SLaurent Vivier 982068f1615SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 983068f1615SLaurent Vivier 984068f1615SLaurent Vivier int32_t compact, n, j, l, m, m1; 985068f1615SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2; 986068f1615SLaurent Vivier 987068f1615SLaurent Vivier aSig = extractFloatx80Frac(a); 988068f1615SLaurent Vivier aExp = extractFloatx80Exp(a); 989068f1615SLaurent Vivier aSign = extractFloatx80Sign(a); 990068f1615SLaurent Vivier 991068f1615SLaurent Vivier if (aExp == 0x7FFF) { 992068f1615SLaurent Vivier if ((uint64_t) (aSig << 1)) { 993068f1615SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 994068f1615SLaurent Vivier } 995068f1615SLaurent Vivier if (aSign) { 996068f1615SLaurent Vivier return packFloatx80(0, 0, 0); 997068f1615SLaurent Vivier } 998068f1615SLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 999068f1615SLaurent Vivier floatx80_infinity.low); 1000068f1615SLaurent Vivier } 1001068f1615SLaurent Vivier 1002068f1615SLaurent Vivier if (aExp == 0 && aSig == 0) { 1003068f1615SLaurent Vivier return packFloatx80(0, one_exp, one_sig); 1004068f1615SLaurent Vivier } 1005068f1615SLaurent Vivier 1006068f1615SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 1007068f1615SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 1008068f1615SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 1009068f1615SLaurent Vivier status->floatx80_rounding_precision = 80; 1010068f1615SLaurent Vivier 1011068f1615SLaurent Vivier fp0 = a; 1012068f1615SLaurent Vivier 1013068f1615SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 1014068f1615SLaurent Vivier 1015068f1615SLaurent Vivier if (compact < 0x3FB98000 || compact > 0x400D80C0) { 1016068f1615SLaurent Vivier /* |X| > 16480 or |X| < 2^(-70) */ 1017068f1615SLaurent Vivier if (compact > 0x3FFF8000) { /* |X| > 16480 */ 1018068f1615SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 1019068f1615SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 1020068f1615SLaurent Vivier 1021068f1615SLaurent Vivier if (aSign) { 1022068f1615SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 1023068f1615SLaurent Vivier 0, -0x1000, aSig, 0, status); 1024068f1615SLaurent Vivier } else { 1025068f1615SLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 1026068f1615SLaurent Vivier 0, 0x8000, aSig, 0, status); 1027068f1615SLaurent Vivier } 1028068f1615SLaurent Vivier } else { /* |X| < 2^(-70) */ 1029068f1615SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 1030068f1615SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 1031068f1615SLaurent Vivier 1032068f1615SLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80( 1033068f1615SLaurent Vivier make_float32(0x3F800000), status), 1034068f1615SLaurent Vivier status); /* 1 + X */ 1035068f1615SLaurent Vivier 1036068f1615SLaurent Vivier float_raise(float_flag_inexact, status); 1037068f1615SLaurent Vivier 1038068f1615SLaurent Vivier return a; 1039068f1615SLaurent Vivier } 1040068f1615SLaurent Vivier } else { /* 2^(-70) <= |X| <= 16480 */ 1041068f1615SLaurent Vivier fp1 = fp0; /* X */ 1042068f1615SLaurent Vivier fp1 = floatx80_mul(fp1, float32_to_floatx80( 1043068f1615SLaurent Vivier make_float32(0x42800000), status), 1044068f1615SLaurent Vivier status); /* X * 64 */ 1045068f1615SLaurent Vivier n = floatx80_to_int32(fp1, status); 1046068f1615SLaurent Vivier fp1 = int32_to_floatx80(n, status); 1047068f1615SLaurent Vivier j = n & 0x3F; 1048068f1615SLaurent Vivier l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */ 1049068f1615SLaurent Vivier if (n < 0 && j) { 1050068f1615SLaurent Vivier /* arithmetic right shift is division and 1051068f1615SLaurent Vivier * round towards minus infinity 1052068f1615SLaurent Vivier */ 1053068f1615SLaurent Vivier l--; 1054068f1615SLaurent Vivier } 1055068f1615SLaurent Vivier m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */ 1056068f1615SLaurent Vivier if (l < 0 && (l & 1)) { 1057068f1615SLaurent Vivier /* arithmetic right shift is division and 1058068f1615SLaurent Vivier * round towards minus infinity 1059068f1615SLaurent Vivier */ 1060068f1615SLaurent Vivier m--; 1061068f1615SLaurent Vivier } 1062068f1615SLaurent Vivier m1 = l - m; 1063068f1615SLaurent Vivier m1 += 0x3FFF; /* ADJFACT IS 2^(M') */ 1064068f1615SLaurent Vivier 1065068f1615SLaurent Vivier adjfact = packFloatx80(0, m1, one_sig); 1066068f1615SLaurent Vivier fact1 = exp2_tbl[j]; 1067068f1615SLaurent Vivier fact1.high += m; 1068068f1615SLaurent Vivier fact2.high = exp2_tbl2[j] >> 16; 1069068f1615SLaurent Vivier fact2.high += m; 1070068f1615SLaurent Vivier fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF); 1071068f1615SLaurent Vivier fact2.low <<= 48; 1072068f1615SLaurent Vivier 1073068f1615SLaurent Vivier fp1 = floatx80_mul(fp1, float32_to_floatx80( 1074068f1615SLaurent Vivier make_float32(0x3C800000), status), 1075068f1615SLaurent Vivier status); /* (1/64)*N */ 1076068f1615SLaurent Vivier fp0 = floatx80_sub(fp0, fp1, status); /* X - (1/64)*INT(64 X) */ 1077068f1615SLaurent Vivier fp2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); /* LOG2 */ 1078068f1615SLaurent Vivier fp0 = floatx80_mul(fp0, fp2, status); /* R */ 1079068f1615SLaurent Vivier 1080068f1615SLaurent Vivier /* EXPR */ 1081068f1615SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 1082068f1615SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2), 1083068f1615SLaurent Vivier status); /* A5 */ 1084068f1615SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C), 1085068f1615SLaurent Vivier status); /* A4 */ 1086068f1615SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */ 1087068f1615SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */ 1088068f1615SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 1089068f1615SLaurent Vivier make_float64(0x3FA5555555554CC1), status), 1090068f1615SLaurent Vivier status); /* A3+S*A5 */ 1091068f1615SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 1092068f1615SLaurent Vivier make_float64(0x3FC5555555554A54), status), 1093068f1615SLaurent Vivier status); /* A2+S*A4 */ 1094068f1615SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */ 1095068f1615SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */ 1096068f1615SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 1097068f1615SLaurent Vivier make_float64(0x3FE0000000000000), status), 1098068f1615SLaurent Vivier status); /* A1+S*(A3+S*A5) */ 1099068f1615SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */ 1100068f1615SLaurent Vivier 1101068f1615SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */ 1102068f1615SLaurent Vivier fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */ 1103068f1615SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */ 1104068f1615SLaurent Vivier 1105068f1615SLaurent Vivier fp0 = floatx80_mul(fp0, fact1, status); 1106068f1615SLaurent Vivier fp0 = floatx80_add(fp0, fact2, status); 1107068f1615SLaurent Vivier fp0 = floatx80_add(fp0, fact1, status); 1108068f1615SLaurent Vivier 1109068f1615SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 1110068f1615SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 1111068f1615SLaurent Vivier 1112068f1615SLaurent Vivier a = floatx80_mul(fp0, adjfact, status); 1113068f1615SLaurent Vivier 1114068f1615SLaurent Vivier float_raise(float_flag_inexact, status); 1115068f1615SLaurent Vivier 1116068f1615SLaurent Vivier return a; 1117068f1615SLaurent Vivier } 1118068f1615SLaurent Vivier } 11196c25be6eSLaurent Vivier 11206c25be6eSLaurent Vivier /*---------------------------------------------------------------------------- 11216c25be6eSLaurent Vivier | 10 to x 11226c25be6eSLaurent Vivier *----------------------------------------------------------------------------*/ 11236c25be6eSLaurent Vivier 11246c25be6eSLaurent Vivier floatx80 floatx80_tentox(floatx80 a, float_status *status) 11256c25be6eSLaurent Vivier { 11266c25be6eSLaurent Vivier flag aSign; 11276c25be6eSLaurent Vivier int32_t aExp; 11286c25be6eSLaurent Vivier uint64_t aSig; 11296c25be6eSLaurent Vivier 11306c25be6eSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 11316c25be6eSLaurent Vivier 11326c25be6eSLaurent Vivier int32_t compact, n, j, l, m, m1; 11336c25be6eSLaurent Vivier floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2; 11346c25be6eSLaurent Vivier 11356c25be6eSLaurent Vivier aSig = extractFloatx80Frac(a); 11366c25be6eSLaurent Vivier aExp = extractFloatx80Exp(a); 11376c25be6eSLaurent Vivier aSign = extractFloatx80Sign(a); 11386c25be6eSLaurent Vivier 11396c25be6eSLaurent Vivier if (aExp == 0x7FFF) { 11406c25be6eSLaurent Vivier if ((uint64_t) (aSig << 1)) { 11416c25be6eSLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 11426c25be6eSLaurent Vivier } 11436c25be6eSLaurent Vivier if (aSign) { 11446c25be6eSLaurent Vivier return packFloatx80(0, 0, 0); 11456c25be6eSLaurent Vivier } 11466c25be6eSLaurent Vivier return packFloatx80(0, floatx80_infinity.high, 11476c25be6eSLaurent Vivier floatx80_infinity.low); 11486c25be6eSLaurent Vivier } 11496c25be6eSLaurent Vivier 11506c25be6eSLaurent Vivier if (aExp == 0 && aSig == 0) { 11516c25be6eSLaurent Vivier return packFloatx80(0, one_exp, one_sig); 11526c25be6eSLaurent Vivier } 11536c25be6eSLaurent Vivier 11546c25be6eSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 11556c25be6eSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 11566c25be6eSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 11576c25be6eSLaurent Vivier status->floatx80_rounding_precision = 80; 11586c25be6eSLaurent Vivier 11596c25be6eSLaurent Vivier fp0 = a; 11606c25be6eSLaurent Vivier 11616c25be6eSLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 11626c25be6eSLaurent Vivier 11636c25be6eSLaurent Vivier if (compact < 0x3FB98000 || compact > 0x400B9B07) { 11646c25be6eSLaurent Vivier /* |X| > 16480 LOG2/LOG10 or |X| < 2^(-70) */ 11656c25be6eSLaurent Vivier if (compact > 0x3FFF8000) { /* |X| > 16480 */ 11666c25be6eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 11676c25be6eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 11686c25be6eSLaurent Vivier 11696c25be6eSLaurent Vivier if (aSign) { 11706c25be6eSLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 11716c25be6eSLaurent Vivier 0, -0x1000, aSig, 0, status); 11726c25be6eSLaurent Vivier } else { 11736c25be6eSLaurent Vivier return roundAndPackFloatx80(status->floatx80_rounding_precision, 11746c25be6eSLaurent Vivier 0, 0x8000, aSig, 0, status); 11756c25be6eSLaurent Vivier } 11766c25be6eSLaurent Vivier } else { /* |X| < 2^(-70) */ 11776c25be6eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 11786c25be6eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 11796c25be6eSLaurent Vivier 11806c25be6eSLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80( 11816c25be6eSLaurent Vivier make_float32(0x3F800000), status), 11826c25be6eSLaurent Vivier status); /* 1 + X */ 11836c25be6eSLaurent Vivier 11846c25be6eSLaurent Vivier float_raise(float_flag_inexact, status); 11856c25be6eSLaurent Vivier 11866c25be6eSLaurent Vivier return a; 11876c25be6eSLaurent Vivier } 11886c25be6eSLaurent Vivier } else { /* 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10 */ 11896c25be6eSLaurent Vivier fp1 = fp0; /* X */ 11906c25be6eSLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 11916c25be6eSLaurent Vivier make_float64(0x406A934F0979A371), 11926c25be6eSLaurent Vivier status), status); /* X*64*LOG10/LOG2 */ 11936c25be6eSLaurent Vivier n = floatx80_to_int32(fp1, status); /* N=INT(X*64*LOG10/LOG2) */ 11946c25be6eSLaurent Vivier fp1 = int32_to_floatx80(n, status); 11956c25be6eSLaurent Vivier 11966c25be6eSLaurent Vivier j = n & 0x3F; 11976c25be6eSLaurent Vivier l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */ 11986c25be6eSLaurent Vivier if (n < 0 && j) { 11996c25be6eSLaurent Vivier /* arithmetic right shift is division and 12006c25be6eSLaurent Vivier * round towards minus infinity 12016c25be6eSLaurent Vivier */ 12026c25be6eSLaurent Vivier l--; 12036c25be6eSLaurent Vivier } 12046c25be6eSLaurent Vivier m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */ 12056c25be6eSLaurent Vivier if (l < 0 && (l & 1)) { 12066c25be6eSLaurent Vivier /* arithmetic right shift is division and 12076c25be6eSLaurent Vivier * round towards minus infinity 12086c25be6eSLaurent Vivier */ 12096c25be6eSLaurent Vivier m--; 12106c25be6eSLaurent Vivier } 12116c25be6eSLaurent Vivier m1 = l - m; 12126c25be6eSLaurent Vivier m1 += 0x3FFF; /* ADJFACT IS 2^(M') */ 12136c25be6eSLaurent Vivier 12146c25be6eSLaurent Vivier adjfact = packFloatx80(0, m1, one_sig); 12156c25be6eSLaurent Vivier fact1 = exp2_tbl[j]; 12166c25be6eSLaurent Vivier fact1.high += m; 12176c25be6eSLaurent Vivier fact2.high = exp2_tbl2[j] >> 16; 12186c25be6eSLaurent Vivier fact2.high += m; 12196c25be6eSLaurent Vivier fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF); 12206c25be6eSLaurent Vivier fact2.low <<= 48; 12216c25be6eSLaurent Vivier 12226c25be6eSLaurent Vivier fp2 = fp1; /* N */ 12236c25be6eSLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 12246c25be6eSLaurent Vivier make_float64(0x3F734413509F8000), status), 12256c25be6eSLaurent Vivier status); /* N*(LOG2/64LOG10)_LEAD */ 12266c25be6eSLaurent Vivier fp3 = packFloatx80(1, 0x3FCD, LIT64(0xC0219DC1DA994FD2)); 12276c25be6eSLaurent Vivier fp2 = floatx80_mul(fp2, fp3, status); /* N*(LOG2/64LOG10)_TRAIL */ 12286c25be6eSLaurent Vivier fp0 = floatx80_sub(fp0, fp1, status); /* X - N L_LEAD */ 12296c25be6eSLaurent Vivier fp0 = floatx80_sub(fp0, fp2, status); /* X - N L_TRAIL */ 12306c25be6eSLaurent Vivier fp2 = packFloatx80(0, 0x4000, LIT64(0x935D8DDDAAA8AC17)); /* LOG10 */ 12316c25be6eSLaurent Vivier fp0 = floatx80_mul(fp0, fp2, status); /* R */ 12326c25be6eSLaurent Vivier 12336c25be6eSLaurent Vivier /* EXPR */ 12346c25be6eSLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 12356c25be6eSLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2), 12366c25be6eSLaurent Vivier status); /* A5 */ 12376c25be6eSLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C), 12386c25be6eSLaurent Vivier status); /* A4 */ 12396c25be6eSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */ 12406c25be6eSLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */ 12416c25be6eSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 12426c25be6eSLaurent Vivier make_float64(0x3FA5555555554CC1), status), 12436c25be6eSLaurent Vivier status); /* A3+S*A5 */ 12446c25be6eSLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 12456c25be6eSLaurent Vivier make_float64(0x3FC5555555554A54), status), 12466c25be6eSLaurent Vivier status); /* A2+S*A4 */ 12476c25be6eSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */ 12486c25be6eSLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */ 12496c25be6eSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 12506c25be6eSLaurent Vivier make_float64(0x3FE0000000000000), status), 12516c25be6eSLaurent Vivier status); /* A1+S*(A3+S*A5) */ 12526c25be6eSLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */ 12536c25be6eSLaurent Vivier 12546c25be6eSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */ 12556c25be6eSLaurent Vivier fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */ 12566c25be6eSLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */ 12576c25be6eSLaurent Vivier 12586c25be6eSLaurent Vivier fp0 = floatx80_mul(fp0, fact1, status); 12596c25be6eSLaurent Vivier fp0 = floatx80_add(fp0, fact2, status); 12606c25be6eSLaurent Vivier fp0 = floatx80_add(fp0, fact1, status); 12616c25be6eSLaurent Vivier 12626c25be6eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 12636c25be6eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 12646c25be6eSLaurent Vivier 12656c25be6eSLaurent Vivier a = floatx80_mul(fp0, adjfact, status); 12666c25be6eSLaurent Vivier 12676c25be6eSLaurent Vivier float_raise(float_flag_inexact, status); 12686c25be6eSLaurent Vivier 12696c25be6eSLaurent Vivier return a; 12706c25be6eSLaurent Vivier } 12716c25be6eSLaurent Vivier } 127227340180SLaurent Vivier 127327340180SLaurent Vivier /*---------------------------------------------------------------------------- 127427340180SLaurent Vivier | Tangent 127527340180SLaurent Vivier *----------------------------------------------------------------------------*/ 127627340180SLaurent Vivier 127727340180SLaurent Vivier floatx80 floatx80_tan(floatx80 a, float_status *status) 127827340180SLaurent Vivier { 127927340180SLaurent Vivier flag aSign, xSign; 128027340180SLaurent Vivier int32_t aExp, xExp; 128127340180SLaurent Vivier uint64_t aSig, xSig; 128227340180SLaurent Vivier 128327340180SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 128427340180SLaurent Vivier 128527340180SLaurent Vivier int32_t compact, l, n, j; 128627340180SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, fp4, fp5, invtwopi, twopi1, twopi2; 128727340180SLaurent Vivier float32 twoto63; 128827340180SLaurent Vivier flag endflag; 128927340180SLaurent Vivier 129027340180SLaurent Vivier aSig = extractFloatx80Frac(a); 129127340180SLaurent Vivier aExp = extractFloatx80Exp(a); 129227340180SLaurent Vivier aSign = extractFloatx80Sign(a); 129327340180SLaurent Vivier 129427340180SLaurent Vivier if (aExp == 0x7FFF) { 129527340180SLaurent Vivier if ((uint64_t) (aSig << 1)) { 129627340180SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 129727340180SLaurent Vivier } 129827340180SLaurent Vivier float_raise(float_flag_invalid, status); 129927340180SLaurent Vivier return floatx80_default_nan(status); 130027340180SLaurent Vivier } 130127340180SLaurent Vivier 130227340180SLaurent Vivier if (aExp == 0 && aSig == 0) { 130327340180SLaurent Vivier return packFloatx80(aSign, 0, 0); 130427340180SLaurent Vivier } 130527340180SLaurent Vivier 130627340180SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 130727340180SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 130827340180SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 130927340180SLaurent Vivier status->floatx80_rounding_precision = 80; 131027340180SLaurent Vivier 131127340180SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 131227340180SLaurent Vivier 131327340180SLaurent Vivier fp0 = a; 131427340180SLaurent Vivier 131527340180SLaurent Vivier if (compact < 0x3FD78000 || compact > 0x4004BC7E) { 131627340180SLaurent Vivier /* 2^(-40) > |X| > 15 PI */ 131727340180SLaurent Vivier if (compact > 0x3FFF8000) { /* |X| >= 15 PI */ 131827340180SLaurent Vivier /* REDUCEX */ 131927340180SLaurent Vivier fp1 = packFloatx80(0, 0, 0); 132027340180SLaurent Vivier if (compact == 0x7FFEFFFF) { 132127340180SLaurent Vivier twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, 132227340180SLaurent Vivier LIT64(0xC90FDAA200000000)); 132327340180SLaurent Vivier twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, 132427340180SLaurent Vivier LIT64(0x85A308D300000000)); 132527340180SLaurent Vivier fp0 = floatx80_add(fp0, twopi1, status); 132627340180SLaurent Vivier fp1 = fp0; 132727340180SLaurent Vivier fp0 = floatx80_add(fp0, twopi2, status); 132827340180SLaurent Vivier fp1 = floatx80_sub(fp1, fp0, status); 132927340180SLaurent Vivier fp1 = floatx80_add(fp1, twopi2, status); 133027340180SLaurent Vivier } 133127340180SLaurent Vivier loop: 133227340180SLaurent Vivier xSign = extractFloatx80Sign(fp0); 133327340180SLaurent Vivier xExp = extractFloatx80Exp(fp0); 133427340180SLaurent Vivier xExp -= 0x3FFF; 133527340180SLaurent Vivier if (xExp <= 28) { 133627340180SLaurent Vivier l = 0; 133727340180SLaurent Vivier endflag = 1; 133827340180SLaurent Vivier } else { 133927340180SLaurent Vivier l = xExp - 27; 134027340180SLaurent Vivier endflag = 0; 134127340180SLaurent Vivier } 134227340180SLaurent Vivier invtwopi = packFloatx80(0, 0x3FFE - l, 134327340180SLaurent Vivier LIT64(0xA2F9836E4E44152A)); /* INVTWOPI */ 134427340180SLaurent Vivier twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); 134527340180SLaurent Vivier twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); 134627340180SLaurent Vivier 134727340180SLaurent Vivier /* SIGN(INARG)*2^63 IN SGL */ 134827340180SLaurent Vivier twoto63 = packFloat32(xSign, 0xBE, 0); 134927340180SLaurent Vivier 135027340180SLaurent Vivier fp2 = floatx80_mul(fp0, invtwopi, status); 135127340180SLaurent Vivier fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), 135227340180SLaurent Vivier status); /* THE FRACT PART OF FP2 IS ROUNDED */ 135327340180SLaurent Vivier fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), 135427340180SLaurent Vivier status); /* FP2 is N */ 135527340180SLaurent Vivier fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */ 135627340180SLaurent Vivier fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */ 135727340180SLaurent Vivier fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */ 135827340180SLaurent Vivier fp4 = floatx80_sub(fp4, fp3, status); /* W-P */ 135927340180SLaurent Vivier fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */ 136027340180SLaurent Vivier fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */ 136127340180SLaurent Vivier fp3 = fp0; /* FP3 is A */ 136227340180SLaurent Vivier fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */ 136327340180SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */ 136427340180SLaurent Vivier 136527340180SLaurent Vivier if (endflag > 0) { 136627340180SLaurent Vivier n = floatx80_to_int32(fp2, status); 136727340180SLaurent Vivier goto tancont; 136827340180SLaurent Vivier } 136927340180SLaurent Vivier fp3 = floatx80_sub(fp3, fp0, status); /* A-R */ 137027340180SLaurent Vivier fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */ 137127340180SLaurent Vivier goto loop; 137227340180SLaurent Vivier } else { 137327340180SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 137427340180SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 137527340180SLaurent Vivier 137627340180SLaurent Vivier a = floatx80_move(a, status); 137727340180SLaurent Vivier 137827340180SLaurent Vivier float_raise(float_flag_inexact, status); 137927340180SLaurent Vivier 138027340180SLaurent Vivier return a; 138127340180SLaurent Vivier } 138227340180SLaurent Vivier } else { 138327340180SLaurent Vivier fp1 = floatx80_mul(fp0, float64_to_floatx80( 138427340180SLaurent Vivier make_float64(0x3FE45F306DC9C883), status), 138527340180SLaurent Vivier status); /* X*2/PI */ 138627340180SLaurent Vivier 138727340180SLaurent Vivier n = floatx80_to_int32(fp1, status); 138827340180SLaurent Vivier j = 32 + n; 138927340180SLaurent Vivier 139027340180SLaurent Vivier fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */ 139127340180SLaurent Vivier fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status), 139227340180SLaurent Vivier status); /* FP0 IS R = (X-Y1)-Y2 */ 139327340180SLaurent Vivier 139427340180SLaurent Vivier tancont: 139527340180SLaurent Vivier if (n & 1) { 139627340180SLaurent Vivier /* NODD */ 139727340180SLaurent Vivier fp1 = fp0; /* R */ 139827340180SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 139927340180SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688), 140027340180SLaurent Vivier status); /* Q4 */ 140127340180SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04), 140227340180SLaurent Vivier status); /* P3 */ 140327340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* SQ4 */ 140427340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* SP3 */ 140527340180SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 140627340180SLaurent Vivier make_float64(0xBF346F59B39BA65F), status), 140727340180SLaurent Vivier status); /* Q3+SQ4 */ 140827340180SLaurent Vivier fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00)); 140927340180SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */ 141027340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* S(Q3+SQ4) */ 141127340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* S(P2+SP3) */ 141227340180SLaurent Vivier fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1)); 141327340180SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */ 141427340180SLaurent Vivier fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA)); 141527340180SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */ 141627340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp0, status); /* S(Q2+S(Q3+SQ4)) */ 141727340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* S(P1+S(P2+SP3)) */ 141827340180SLaurent Vivier fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE)); 141927340180SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */ 142027340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* RS(P1+S(P2+SP3)) */ 142127340180SLaurent Vivier fp0 = floatx80_mul(fp0, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */ 142227340180SLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); /* R+RS(P1+S(P2+SP3)) */ 142327340180SLaurent Vivier fp0 = floatx80_add(fp0, float32_to_floatx80( 142427340180SLaurent Vivier make_float32(0x3F800000), status), 142527340180SLaurent Vivier status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */ 142627340180SLaurent Vivier 142727340180SLaurent Vivier xSign = extractFloatx80Sign(fp1); 142827340180SLaurent Vivier xExp = extractFloatx80Exp(fp1); 142927340180SLaurent Vivier xSig = extractFloatx80Frac(fp1); 143027340180SLaurent Vivier xSign ^= 1; 143127340180SLaurent Vivier fp1 = packFloatx80(xSign, xExp, xSig); 143227340180SLaurent Vivier 143327340180SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 143427340180SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 143527340180SLaurent Vivier 143627340180SLaurent Vivier a = floatx80_div(fp0, fp1, status); 143727340180SLaurent Vivier 143827340180SLaurent Vivier float_raise(float_flag_inexact, status); 143927340180SLaurent Vivier 144027340180SLaurent Vivier return a; 144127340180SLaurent Vivier } else { 144227340180SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */ 144327340180SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688), 144427340180SLaurent Vivier status); /* Q4 */ 144527340180SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04), 144627340180SLaurent Vivier status); /* P3 */ 144727340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* SQ4 */ 144827340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* SP3 */ 144927340180SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 145027340180SLaurent Vivier make_float64(0xBF346F59B39BA65F), status), 145127340180SLaurent Vivier status); /* Q3+SQ4 */ 145227340180SLaurent Vivier fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00)); 145327340180SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */ 145427340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S(Q3+SQ4) */ 145527340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S(P2+SP3) */ 145627340180SLaurent Vivier fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1)); 145727340180SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */ 145827340180SLaurent Vivier fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA)); 145927340180SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */ 146027340180SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* S(Q2+S(Q3+SQ4)) */ 146127340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* S(P1+S(P2+SP3)) */ 146227340180SLaurent Vivier fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE)); 146327340180SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */ 146427340180SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* RS(P1+S(P2+SP3)) */ 146527340180SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */ 146627340180SLaurent Vivier fp0 = floatx80_add(fp0, fp2, status); /* R+RS(P1+S(P2+SP3)) */ 146727340180SLaurent Vivier fp1 = floatx80_add(fp1, float32_to_floatx80( 146827340180SLaurent Vivier make_float32(0x3F800000), status), 146927340180SLaurent Vivier status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */ 147027340180SLaurent Vivier 147127340180SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 147227340180SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 147327340180SLaurent Vivier 147427340180SLaurent Vivier a = floatx80_div(fp0, fp1, status); 147527340180SLaurent Vivier 147627340180SLaurent Vivier float_raise(float_flag_inexact, status); 147727340180SLaurent Vivier 147827340180SLaurent Vivier return a; 147927340180SLaurent Vivier } 148027340180SLaurent Vivier } 148127340180SLaurent Vivier } 14825add1ac4SLaurent Vivier 14835add1ac4SLaurent Vivier /*---------------------------------------------------------------------------- 14845add1ac4SLaurent Vivier | Sine 14855add1ac4SLaurent Vivier *----------------------------------------------------------------------------*/ 14865add1ac4SLaurent Vivier 14875add1ac4SLaurent Vivier floatx80 floatx80_sin(floatx80 a, float_status *status) 14885add1ac4SLaurent Vivier { 14895add1ac4SLaurent Vivier flag aSign, xSign; 14905add1ac4SLaurent Vivier int32_t aExp, xExp; 14915add1ac4SLaurent Vivier uint64_t aSig, xSig; 14925add1ac4SLaurent Vivier 14935add1ac4SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 14945add1ac4SLaurent Vivier 14955add1ac4SLaurent Vivier int32_t compact, l, n, j; 14965add1ac4SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2; 14975add1ac4SLaurent Vivier float32 posneg1, twoto63; 14985add1ac4SLaurent Vivier flag adjn, endflag; 14995add1ac4SLaurent Vivier 15005add1ac4SLaurent Vivier aSig = extractFloatx80Frac(a); 15015add1ac4SLaurent Vivier aExp = extractFloatx80Exp(a); 15025add1ac4SLaurent Vivier aSign = extractFloatx80Sign(a); 15035add1ac4SLaurent Vivier 15045add1ac4SLaurent Vivier if (aExp == 0x7FFF) { 15055add1ac4SLaurent Vivier if ((uint64_t) (aSig << 1)) { 15065add1ac4SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 15075add1ac4SLaurent Vivier } 15085add1ac4SLaurent Vivier float_raise(float_flag_invalid, status); 15095add1ac4SLaurent Vivier return floatx80_default_nan(status); 15105add1ac4SLaurent Vivier } 15115add1ac4SLaurent Vivier 15125add1ac4SLaurent Vivier if (aExp == 0 && aSig == 0) { 15135add1ac4SLaurent Vivier return packFloatx80(aSign, 0, 0); 15145add1ac4SLaurent Vivier } 15155add1ac4SLaurent Vivier 15165add1ac4SLaurent Vivier adjn = 0; 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 if (adjn) { 15935add1ac4SLaurent Vivier /* COSTINY */ 15945add1ac4SLaurent Vivier a = floatx80_sub(fp0, float32_to_floatx80( 15955add1ac4SLaurent Vivier make_float32(0x00800000), status), status); 15965add1ac4SLaurent Vivier } else { 15975add1ac4SLaurent Vivier /* SINTINY */ 15985add1ac4SLaurent Vivier a = floatx80_move(a, status); 15995add1ac4SLaurent Vivier } 16005add1ac4SLaurent Vivier float_raise(float_flag_inexact, status); 16015add1ac4SLaurent Vivier 16025add1ac4SLaurent Vivier return a; 16035add1ac4SLaurent Vivier } 16045add1ac4SLaurent Vivier } else { 16055add1ac4SLaurent Vivier fp1 = floatx80_mul(fp0, float64_to_floatx80( 16065add1ac4SLaurent Vivier make_float64(0x3FE45F306DC9C883), status), 16075add1ac4SLaurent Vivier status); /* X*2/PI */ 16085add1ac4SLaurent Vivier 16095add1ac4SLaurent Vivier n = floatx80_to_int32(fp1, status); 16105add1ac4SLaurent Vivier j = 32 + n; 16115add1ac4SLaurent Vivier 16125add1ac4SLaurent Vivier fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */ 16135add1ac4SLaurent Vivier fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status), 16145add1ac4SLaurent Vivier status); /* FP0 IS R = (X-Y1)-Y2 */ 16155add1ac4SLaurent Vivier 16165add1ac4SLaurent Vivier sincont: 16175add1ac4SLaurent Vivier if ((n + adjn) & 1) { 16185add1ac4SLaurent Vivier /* COSPOLY */ 16195add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */ 16205add1ac4SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */ 16215add1ac4SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3), 16225add1ac4SLaurent Vivier status); /* B8 */ 16235add1ac4SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19), 16245add1ac4SLaurent Vivier status); /* B7 */ 16255add1ac4SLaurent Vivier 16265add1ac4SLaurent Vivier xSign = extractFloatx80Sign(fp0); /* X IS S */ 16275add1ac4SLaurent Vivier xExp = extractFloatx80Exp(fp0); 16285add1ac4SLaurent Vivier xSig = extractFloatx80Frac(fp0); 16295add1ac4SLaurent Vivier 16305add1ac4SLaurent Vivier if (((n + adjn) >> 1) & 1) { 16315add1ac4SLaurent Vivier xSign ^= 1; 16325add1ac4SLaurent Vivier posneg1 = make_float32(0xBF800000); /* -1 */ 16335add1ac4SLaurent Vivier } else { 16345add1ac4SLaurent Vivier xSign ^= 0; 16355add1ac4SLaurent Vivier posneg1 = make_float32(0x3F800000); /* 1 */ 16365add1ac4SLaurent Vivier } /* X IS NOW R'= SGN*R */ 16375add1ac4SLaurent Vivier 16385add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */ 16395add1ac4SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */ 16405add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 16415add1ac4SLaurent Vivier make_float64(0x3E21EED90612C972), status), 16425add1ac4SLaurent Vivier status); /* B6+TB8 */ 16435add1ac4SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 16445add1ac4SLaurent Vivier make_float64(0xBE927E4FB79D9FCF), status), 16455add1ac4SLaurent Vivier status); /* B5+TB7 */ 16465add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */ 16475add1ac4SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */ 16485add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 16495add1ac4SLaurent Vivier make_float64(0x3EFA01A01A01D423), status), 16505add1ac4SLaurent Vivier status); /* B4+T(B6+TB8) */ 16515add1ac4SLaurent Vivier fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438)); 16525add1ac4SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */ 16535add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */ 16545add1ac4SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */ 16555add1ac4SLaurent Vivier fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E)); 16565add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */ 16575add1ac4SLaurent Vivier fp1 = floatx80_add(fp1, float32_to_floatx80( 16585add1ac4SLaurent Vivier make_float32(0xBF000000), status), 16595add1ac4SLaurent Vivier status); /* B1+T(B3+T(B5+TB7)) */ 16605add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */ 16615add1ac4SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* [B1+T(B3+T(B5+TB7))]+ 16625add1ac4SLaurent Vivier * [S(B2+T(B4+T(B6+TB8)))] 16635add1ac4SLaurent Vivier */ 16645add1ac4SLaurent Vivier 16655add1ac4SLaurent Vivier x = packFloatx80(xSign, xExp, xSig); 16665add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, x, status); 16675add1ac4SLaurent Vivier 16685add1ac4SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 16695add1ac4SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 16705add1ac4SLaurent Vivier 16715add1ac4SLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status); 16725add1ac4SLaurent Vivier 16735add1ac4SLaurent Vivier float_raise(float_flag_inexact, status); 16745add1ac4SLaurent Vivier 16755add1ac4SLaurent Vivier return a; 16765add1ac4SLaurent Vivier } else { 16775add1ac4SLaurent Vivier /* SINPOLY */ 16785add1ac4SLaurent Vivier xSign = extractFloatx80Sign(fp0); /* X IS R */ 16795add1ac4SLaurent Vivier xExp = extractFloatx80Exp(fp0); 16805add1ac4SLaurent Vivier xSig = extractFloatx80Frac(fp0); 16815add1ac4SLaurent Vivier 16825add1ac4SLaurent Vivier xSign ^= ((n + adjn) >> 1) & 1; /* X IS NOW R'= SGN*R */ 16835add1ac4SLaurent Vivier 16845add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */ 16855add1ac4SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */ 16865add1ac4SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5), 16875add1ac4SLaurent Vivier status); /* A7 */ 16885add1ac4SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1), 16895add1ac4SLaurent Vivier status); /* A6 */ 16905add1ac4SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */ 16915add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */ 16925add1ac4SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 16935add1ac4SLaurent Vivier make_float64(0xBE5AE6452A118AE4), status), 16945add1ac4SLaurent Vivier status); /* A5+T*A7 */ 16955add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 16965add1ac4SLaurent Vivier make_float64(0x3EC71DE3A5341531), status), 16975add1ac4SLaurent Vivier status); /* A4+T*A6 */ 16985add1ac4SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */ 16995add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */ 17005add1ac4SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 17015add1ac4SLaurent Vivier make_float64(0xBF2A01A01A018B59), status), 17025add1ac4SLaurent Vivier status); /* A3+T(A5+TA7) */ 17035add1ac4SLaurent Vivier fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF)); 17045add1ac4SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */ 17055add1ac4SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */ 17065add1ac4SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */ 17075add1ac4SLaurent Vivier fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99)); 17085add1ac4SLaurent Vivier fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */ 17095add1ac4SLaurent Vivier fp1 = floatx80_add(fp1, fp2, 17105add1ac4SLaurent Vivier status); /* [A1+T(A3+T(A5+TA7))]+ 17115add1ac4SLaurent Vivier * [S(A2+T(A4+TA6))] 17125add1ac4SLaurent Vivier */ 17135add1ac4SLaurent Vivier 17145add1ac4SLaurent Vivier x = packFloatx80(xSign, xExp, xSig); 17155add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, x, status); /* R'*S */ 17165add1ac4SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */ 17175add1ac4SLaurent Vivier 17185add1ac4SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 17195add1ac4SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 17205add1ac4SLaurent Vivier 17215add1ac4SLaurent Vivier a = floatx80_add(fp0, x, status); 17225add1ac4SLaurent Vivier 17235add1ac4SLaurent Vivier float_raise(float_flag_inexact, status); 17245add1ac4SLaurent Vivier 17255add1ac4SLaurent Vivier return a; 17265add1ac4SLaurent Vivier } 17275add1ac4SLaurent Vivier } 17285add1ac4SLaurent Vivier } 172968d0ed37SLaurent Vivier 173068d0ed37SLaurent Vivier /*---------------------------------------------------------------------------- 173168d0ed37SLaurent Vivier | Cosine 173268d0ed37SLaurent Vivier *----------------------------------------------------------------------------*/ 173368d0ed37SLaurent Vivier 173468d0ed37SLaurent Vivier floatx80 floatx80_cos(floatx80 a, float_status *status) 173568d0ed37SLaurent Vivier { 173668d0ed37SLaurent Vivier flag aSign, xSign; 173768d0ed37SLaurent Vivier int32_t aExp, xExp; 173868d0ed37SLaurent Vivier uint64_t aSig, xSig; 173968d0ed37SLaurent Vivier 174068d0ed37SLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 174168d0ed37SLaurent Vivier 174268d0ed37SLaurent Vivier int32_t compact, l, n, j; 174368d0ed37SLaurent Vivier floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2; 174468d0ed37SLaurent Vivier float32 posneg1, twoto63; 174568d0ed37SLaurent Vivier flag adjn, endflag; 174668d0ed37SLaurent Vivier 174768d0ed37SLaurent Vivier aSig = extractFloatx80Frac(a); 174868d0ed37SLaurent Vivier aExp = extractFloatx80Exp(a); 174968d0ed37SLaurent Vivier aSign = extractFloatx80Sign(a); 175068d0ed37SLaurent Vivier 175168d0ed37SLaurent Vivier if (aExp == 0x7FFF) { 175268d0ed37SLaurent Vivier if ((uint64_t) (aSig << 1)) { 175368d0ed37SLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 175468d0ed37SLaurent Vivier } 175568d0ed37SLaurent Vivier float_raise(float_flag_invalid, status); 175668d0ed37SLaurent Vivier return floatx80_default_nan(status); 175768d0ed37SLaurent Vivier } 175868d0ed37SLaurent Vivier 175968d0ed37SLaurent Vivier if (aExp == 0 && aSig == 0) { 176068d0ed37SLaurent Vivier return packFloatx80(0, one_exp, one_sig); 176168d0ed37SLaurent Vivier } 176268d0ed37SLaurent Vivier 176368d0ed37SLaurent Vivier adjn = 1; 176468d0ed37SLaurent Vivier 176568d0ed37SLaurent Vivier user_rnd_mode = status->float_rounding_mode; 176668d0ed37SLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 176768d0ed37SLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 176868d0ed37SLaurent Vivier status->floatx80_rounding_precision = 80; 176968d0ed37SLaurent Vivier 177068d0ed37SLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 177168d0ed37SLaurent Vivier 177268d0ed37SLaurent Vivier fp0 = a; 177368d0ed37SLaurent Vivier 177468d0ed37SLaurent Vivier if (compact < 0x3FD78000 || compact > 0x4004BC7E) { 177568d0ed37SLaurent Vivier /* 2^(-40) > |X| > 15 PI */ 177668d0ed37SLaurent Vivier if (compact > 0x3FFF8000) { /* |X| >= 15 PI */ 177768d0ed37SLaurent Vivier /* REDUCEX */ 177868d0ed37SLaurent Vivier fp1 = packFloatx80(0, 0, 0); 177968d0ed37SLaurent Vivier if (compact == 0x7FFEFFFF) { 178068d0ed37SLaurent Vivier twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, 178168d0ed37SLaurent Vivier LIT64(0xC90FDAA200000000)); 178268d0ed37SLaurent Vivier twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, 178368d0ed37SLaurent Vivier LIT64(0x85A308D300000000)); 178468d0ed37SLaurent Vivier fp0 = floatx80_add(fp0, twopi1, status); 178568d0ed37SLaurent Vivier fp1 = fp0; 178668d0ed37SLaurent Vivier fp0 = floatx80_add(fp0, twopi2, status); 178768d0ed37SLaurent Vivier fp1 = floatx80_sub(fp1, fp0, status); 178868d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, twopi2, status); 178968d0ed37SLaurent Vivier } 179068d0ed37SLaurent Vivier loop: 179168d0ed37SLaurent Vivier xSign = extractFloatx80Sign(fp0); 179268d0ed37SLaurent Vivier xExp = extractFloatx80Exp(fp0); 179368d0ed37SLaurent Vivier xExp -= 0x3FFF; 179468d0ed37SLaurent Vivier if (xExp <= 28) { 179568d0ed37SLaurent Vivier l = 0; 179668d0ed37SLaurent Vivier endflag = 1; 179768d0ed37SLaurent Vivier } else { 179868d0ed37SLaurent Vivier l = xExp - 27; 179968d0ed37SLaurent Vivier endflag = 0; 180068d0ed37SLaurent Vivier } 180168d0ed37SLaurent Vivier invtwopi = packFloatx80(0, 0x3FFE - l, 180268d0ed37SLaurent Vivier LIT64(0xA2F9836E4E44152A)); /* INVTWOPI */ 180368d0ed37SLaurent Vivier twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); 180468d0ed37SLaurent Vivier twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); 180568d0ed37SLaurent Vivier 180668d0ed37SLaurent Vivier /* SIGN(INARG)*2^63 IN SGL */ 180768d0ed37SLaurent Vivier twoto63 = packFloat32(xSign, 0xBE, 0); 180868d0ed37SLaurent Vivier 180968d0ed37SLaurent Vivier fp2 = floatx80_mul(fp0, invtwopi, status); 181068d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), 181168d0ed37SLaurent Vivier status); /* THE FRACT PART OF FP2 IS ROUNDED */ 181268d0ed37SLaurent Vivier fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), 181368d0ed37SLaurent Vivier status); /* FP2 is N */ 181468d0ed37SLaurent Vivier fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */ 181568d0ed37SLaurent Vivier fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */ 181668d0ed37SLaurent Vivier fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */ 181768d0ed37SLaurent Vivier fp4 = floatx80_sub(fp4, fp3, status); /* W-P */ 181868d0ed37SLaurent Vivier fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */ 181968d0ed37SLaurent Vivier fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */ 182068d0ed37SLaurent Vivier fp3 = fp0; /* FP3 is A */ 182168d0ed37SLaurent Vivier fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */ 182268d0ed37SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */ 182368d0ed37SLaurent Vivier 182468d0ed37SLaurent Vivier if (endflag > 0) { 182568d0ed37SLaurent Vivier n = floatx80_to_int32(fp2, status); 182668d0ed37SLaurent Vivier goto sincont; 182768d0ed37SLaurent Vivier } 182868d0ed37SLaurent Vivier fp3 = floatx80_sub(fp3, fp0, status); /* A-R */ 182968d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */ 183068d0ed37SLaurent Vivier goto loop; 183168d0ed37SLaurent Vivier } else { 183268d0ed37SLaurent Vivier /* SINSM */ 183368d0ed37SLaurent Vivier fp0 = float32_to_floatx80(make_float32(0x3F800000), status); /* 1 */ 183468d0ed37SLaurent Vivier 183568d0ed37SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 183668d0ed37SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 183768d0ed37SLaurent Vivier 183868d0ed37SLaurent Vivier if (adjn) { 183968d0ed37SLaurent Vivier /* COSTINY */ 184068d0ed37SLaurent Vivier a = floatx80_sub(fp0, float32_to_floatx80( 184168d0ed37SLaurent Vivier make_float32(0x00800000), status), 184268d0ed37SLaurent Vivier status); 184368d0ed37SLaurent Vivier } else { 184468d0ed37SLaurent Vivier /* SINTINY */ 184568d0ed37SLaurent Vivier a = floatx80_move(a, status); 184668d0ed37SLaurent Vivier } 184768d0ed37SLaurent Vivier float_raise(float_flag_inexact, status); 184868d0ed37SLaurent Vivier 184968d0ed37SLaurent Vivier return a; 185068d0ed37SLaurent Vivier } 185168d0ed37SLaurent Vivier } else { 185268d0ed37SLaurent Vivier fp1 = floatx80_mul(fp0, float64_to_floatx80( 185368d0ed37SLaurent Vivier make_float64(0x3FE45F306DC9C883), status), 185468d0ed37SLaurent Vivier status); /* X*2/PI */ 185568d0ed37SLaurent Vivier 185668d0ed37SLaurent Vivier n = floatx80_to_int32(fp1, status); 185768d0ed37SLaurent Vivier j = 32 + n; 185868d0ed37SLaurent Vivier 185968d0ed37SLaurent Vivier fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */ 186068d0ed37SLaurent Vivier fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status), 186168d0ed37SLaurent Vivier status); /* FP0 IS R = (X-Y1)-Y2 */ 186268d0ed37SLaurent Vivier 186368d0ed37SLaurent Vivier sincont: 186468d0ed37SLaurent Vivier if ((n + adjn) & 1) { 186568d0ed37SLaurent Vivier /* COSPOLY */ 186668d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */ 186768d0ed37SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */ 186868d0ed37SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3), 186968d0ed37SLaurent Vivier status); /* B8 */ 187068d0ed37SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19), 187168d0ed37SLaurent Vivier status); /* B7 */ 187268d0ed37SLaurent Vivier 187368d0ed37SLaurent Vivier xSign = extractFloatx80Sign(fp0); /* X IS S */ 187468d0ed37SLaurent Vivier xExp = extractFloatx80Exp(fp0); 187568d0ed37SLaurent Vivier xSig = extractFloatx80Frac(fp0); 187668d0ed37SLaurent Vivier 187768d0ed37SLaurent Vivier if (((n + adjn) >> 1) & 1) { 187868d0ed37SLaurent Vivier xSign ^= 1; 187968d0ed37SLaurent Vivier posneg1 = make_float32(0xBF800000); /* -1 */ 188068d0ed37SLaurent Vivier } else { 188168d0ed37SLaurent Vivier xSign ^= 0; 188268d0ed37SLaurent Vivier posneg1 = make_float32(0x3F800000); /* 1 */ 188368d0ed37SLaurent Vivier } /* X IS NOW R'= SGN*R */ 188468d0ed37SLaurent Vivier 188568d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */ 188668d0ed37SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */ 188768d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 188868d0ed37SLaurent Vivier make_float64(0x3E21EED90612C972), status), 188968d0ed37SLaurent Vivier status); /* B6+TB8 */ 189068d0ed37SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 189168d0ed37SLaurent Vivier make_float64(0xBE927E4FB79D9FCF), status), 189268d0ed37SLaurent Vivier status); /* B5+TB7 */ 189368d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */ 189468d0ed37SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */ 189568d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 189668d0ed37SLaurent Vivier make_float64(0x3EFA01A01A01D423), status), 189768d0ed37SLaurent Vivier status); /* B4+T(B6+TB8) */ 189868d0ed37SLaurent Vivier fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438)); 189968d0ed37SLaurent Vivier fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */ 190068d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */ 190168d0ed37SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */ 190268d0ed37SLaurent Vivier fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E)); 190368d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */ 190468d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, float32_to_floatx80( 190568d0ed37SLaurent Vivier make_float32(0xBF000000), status), 190668d0ed37SLaurent Vivier status); /* B1+T(B3+T(B5+TB7)) */ 190768d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */ 190868d0ed37SLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); 190968d0ed37SLaurent Vivier /* [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))] */ 191068d0ed37SLaurent Vivier 191168d0ed37SLaurent Vivier x = packFloatx80(xSign, xExp, xSig); 191268d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, x, status); 191368d0ed37SLaurent Vivier 191468d0ed37SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 191568d0ed37SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 191668d0ed37SLaurent Vivier 191768d0ed37SLaurent Vivier a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status); 191868d0ed37SLaurent Vivier 191968d0ed37SLaurent Vivier float_raise(float_flag_inexact, status); 192068d0ed37SLaurent Vivier 192168d0ed37SLaurent Vivier return a; 192268d0ed37SLaurent Vivier } else { 192368d0ed37SLaurent Vivier /* SINPOLY */ 192468d0ed37SLaurent Vivier xSign = extractFloatx80Sign(fp0); /* X IS R */ 192568d0ed37SLaurent Vivier xExp = extractFloatx80Exp(fp0); 192668d0ed37SLaurent Vivier xSig = extractFloatx80Frac(fp0); 192768d0ed37SLaurent Vivier 192868d0ed37SLaurent Vivier xSign ^= ((n + adjn) >> 1) & 1; /* X IS NOW R'= SGN*R */ 192968d0ed37SLaurent Vivier 193068d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */ 193168d0ed37SLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */ 193268d0ed37SLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5), 193368d0ed37SLaurent Vivier status); /* A7 */ 193468d0ed37SLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1), 193568d0ed37SLaurent Vivier status); /* A6 */ 193668d0ed37SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */ 193768d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */ 193868d0ed37SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 193968d0ed37SLaurent Vivier make_float64(0xBE5AE6452A118AE4), status), 194068d0ed37SLaurent Vivier status); /* A5+T*A7 */ 194168d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 194268d0ed37SLaurent Vivier make_float64(0x3EC71DE3A5341531), status), 194368d0ed37SLaurent Vivier status); /* A4+T*A6 */ 194468d0ed37SLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */ 194568d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */ 194668d0ed37SLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 194768d0ed37SLaurent Vivier make_float64(0xBF2A01A01A018B59), status), 194868d0ed37SLaurent Vivier status); /* A3+T(A5+TA7) */ 194968d0ed37SLaurent Vivier fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF)); 195068d0ed37SLaurent Vivier fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */ 195168d0ed37SLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */ 195268d0ed37SLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */ 195368d0ed37SLaurent Vivier fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99)); 195468d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */ 195568d0ed37SLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); 195668d0ed37SLaurent Vivier /* [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))] */ 195768d0ed37SLaurent Vivier 195868d0ed37SLaurent Vivier x = packFloatx80(xSign, xExp, xSig); 195968d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, x, status); /* R'*S */ 196068d0ed37SLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */ 196168d0ed37SLaurent Vivier 196268d0ed37SLaurent Vivier status->float_rounding_mode = user_rnd_mode; 196368d0ed37SLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 196468d0ed37SLaurent Vivier 196568d0ed37SLaurent Vivier a = floatx80_add(fp0, x, status); 196668d0ed37SLaurent Vivier 196768d0ed37SLaurent Vivier float_raise(float_flag_inexact, status); 196868d0ed37SLaurent Vivier 196968d0ed37SLaurent Vivier return a; 197068d0ed37SLaurent Vivier } 197168d0ed37SLaurent Vivier } 197268d0ed37SLaurent Vivier } 19738c992abcSLaurent Vivier 19748c992abcSLaurent Vivier /*---------------------------------------------------------------------------- 19758c992abcSLaurent Vivier | Arc tangent 19768c992abcSLaurent Vivier *----------------------------------------------------------------------------*/ 19778c992abcSLaurent Vivier 19788c992abcSLaurent Vivier floatx80 floatx80_atan(floatx80 a, float_status *status) 19798c992abcSLaurent Vivier { 19808c992abcSLaurent Vivier flag aSign; 19818c992abcSLaurent Vivier int32_t aExp; 19828c992abcSLaurent Vivier uint64_t aSig; 19838c992abcSLaurent Vivier 19848c992abcSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 19858c992abcSLaurent Vivier 19868c992abcSLaurent Vivier int32_t compact, tbl_index; 19878c992abcSLaurent Vivier floatx80 fp0, fp1, fp2, fp3, xsave; 19888c992abcSLaurent Vivier 19898c992abcSLaurent Vivier aSig = extractFloatx80Frac(a); 19908c992abcSLaurent Vivier aExp = extractFloatx80Exp(a); 19918c992abcSLaurent Vivier aSign = extractFloatx80Sign(a); 19928c992abcSLaurent Vivier 19938c992abcSLaurent Vivier if (aExp == 0x7FFF) { 19948c992abcSLaurent Vivier if ((uint64_t) (aSig << 1)) { 19958c992abcSLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 19968c992abcSLaurent Vivier } 19978c992abcSLaurent Vivier a = packFloatx80(aSign, piby2_exp, pi_sig); 19988c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 19998c992abcSLaurent Vivier return floatx80_move(a, status); 20008c992abcSLaurent Vivier } 20018c992abcSLaurent Vivier 20028c992abcSLaurent Vivier if (aExp == 0 && aSig == 0) { 20038c992abcSLaurent Vivier return packFloatx80(aSign, 0, 0); 20048c992abcSLaurent Vivier } 20058c992abcSLaurent Vivier 20068c992abcSLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 20078c992abcSLaurent Vivier 20088c992abcSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 20098c992abcSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 20108c992abcSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 20118c992abcSLaurent Vivier status->floatx80_rounding_precision = 80; 20128c992abcSLaurent Vivier 20138c992abcSLaurent Vivier if (compact < 0x3FFB8000 || compact > 0x4002FFFF) { 20148c992abcSLaurent Vivier /* |X| >= 16 or |X| < 1/16 */ 20158c992abcSLaurent Vivier if (compact > 0x3FFF8000) { /* |X| >= 16 */ 20168c992abcSLaurent Vivier if (compact > 0x40638000) { /* |X| > 2^(100) */ 20178c992abcSLaurent Vivier fp0 = packFloatx80(aSign, piby2_exp, pi_sig); 20188c992abcSLaurent Vivier fp1 = packFloatx80(aSign, 0x0001, one_sig); 20198c992abcSLaurent Vivier 20208c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 20218c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 20228c992abcSLaurent Vivier 20238c992abcSLaurent Vivier a = floatx80_sub(fp0, fp1, status); 20248c992abcSLaurent Vivier 20258c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 20268c992abcSLaurent Vivier 20278c992abcSLaurent Vivier return a; 20288c992abcSLaurent Vivier } else { 20298c992abcSLaurent Vivier fp0 = a; 20308c992abcSLaurent Vivier fp1 = packFloatx80(1, one_exp, one_sig); /* -1 */ 20318c992abcSLaurent Vivier fp1 = floatx80_div(fp1, fp0, status); /* X' = -1/X */ 20328c992abcSLaurent Vivier xsave = fp1; 20338c992abcSLaurent Vivier fp0 = floatx80_mul(fp1, fp1, status); /* Y = X'*X' */ 20348c992abcSLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */ 20358c992abcSLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBFB70BF398539E6A), 20368c992abcSLaurent Vivier status); /* C5 */ 20378c992abcSLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3FBC7187962D1D7D), 20388c992abcSLaurent Vivier status); /* C4 */ 20398c992abcSLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* Z*C5 */ 20408c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* Z*C4 */ 20418c992abcSLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 20428c992abcSLaurent Vivier make_float64(0xBFC24924827107B8), status), 20438c992abcSLaurent Vivier status); /* C3+Z*C5 */ 20448c992abcSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 20458c992abcSLaurent Vivier make_float64(0x3FC999999996263E), status), 20468c992abcSLaurent Vivier status); /* C2+Z*C4 */ 20478c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* Z*(C3+Z*C5) */ 20488c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* Y*(C2+Z*C4) */ 20498c992abcSLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 20508c992abcSLaurent Vivier make_float64(0xBFD5555555555536), status), 20518c992abcSLaurent Vivier status); /* C1+Z*(C3+Z*C5) */ 20528c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, xsave, status); /* X'*Y */ 20538c992abcSLaurent Vivier /* [Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)] */ 20548c992abcSLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); 20558c992abcSLaurent Vivier /* X'*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) ?? */ 20568c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); 20578c992abcSLaurent Vivier fp0 = floatx80_add(fp0, xsave, status); 20588c992abcSLaurent Vivier fp1 = packFloatx80(aSign, piby2_exp, pi_sig); 20598c992abcSLaurent Vivier 20608c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 20618c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 20628c992abcSLaurent Vivier 20638c992abcSLaurent Vivier a = floatx80_add(fp0, fp1, status); 20648c992abcSLaurent Vivier 20658c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 20668c992abcSLaurent Vivier 20678c992abcSLaurent Vivier return a; 20688c992abcSLaurent Vivier } 20698c992abcSLaurent Vivier } else { /* |X| < 1/16 */ 20708c992abcSLaurent Vivier if (compact < 0x3FD78000) { /* |X| < 2^(-40) */ 20718c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 20728c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 20738c992abcSLaurent Vivier 20748c992abcSLaurent Vivier a = floatx80_move(a, status); 20758c992abcSLaurent Vivier 20768c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 20778c992abcSLaurent Vivier 20788c992abcSLaurent Vivier return a; 20798c992abcSLaurent Vivier } else { 20808c992abcSLaurent Vivier fp0 = a; 20818c992abcSLaurent Vivier xsave = a; 20828c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, fp0, status); /* Y = X*X */ 20838c992abcSLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */ 20848c992abcSLaurent Vivier fp2 = float64_to_floatx80(make_float64(0x3FB344447F876989), 20858c992abcSLaurent Vivier status); /* B6 */ 20868c992abcSLaurent Vivier fp3 = float64_to_floatx80(make_float64(0xBFB744EE7FAF45DB), 20878c992abcSLaurent Vivier status); /* B5 */ 20888c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* Z*B6 */ 20898c992abcSLaurent Vivier fp3 = floatx80_mul(fp3, fp1, status); /* Z*B5 */ 20908c992abcSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 20918c992abcSLaurent Vivier make_float64(0x3FBC71C646940220), status), 20928c992abcSLaurent Vivier status); /* B4+Z*B6 */ 20938c992abcSLaurent Vivier fp3 = floatx80_add(fp3, float64_to_floatx80( 20948c992abcSLaurent Vivier make_float64(0xBFC24924921872F9), 20958c992abcSLaurent Vivier status), status); /* B3+Z*B5 */ 20968c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* Z*(B4+Z*B6) */ 20978c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, fp3, status); /* Z*(B3+Z*B5) */ 20988c992abcSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 20998c992abcSLaurent Vivier make_float64(0x3FC9999999998FA9), status), 21008c992abcSLaurent Vivier status); /* B2+Z*(B4+Z*B6) */ 21018c992abcSLaurent Vivier fp1 = floatx80_add(fp1, float64_to_floatx80( 21028c992abcSLaurent Vivier make_float64(0xBFD5555555555555), status), 21038c992abcSLaurent Vivier status); /* B1+Z*(B3+Z*B5) */ 21048c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp0, status); /* Y*(B2+Z*(B4+Z*B6)) */ 21058c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, xsave, status); /* X*Y */ 21068c992abcSLaurent Vivier /* [B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))] */ 21078c992abcSLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); 21088c992abcSLaurent Vivier /* X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) */ 21098c992abcSLaurent Vivier fp0 = floatx80_mul(fp0, fp1, status); 21108c992abcSLaurent Vivier 21118c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 21128c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 21138c992abcSLaurent Vivier 21148c992abcSLaurent Vivier a = floatx80_add(fp0, xsave, status); 21158c992abcSLaurent Vivier 21168c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 21178c992abcSLaurent Vivier 21188c992abcSLaurent Vivier return a; 21198c992abcSLaurent Vivier } 21208c992abcSLaurent Vivier } 21218c992abcSLaurent Vivier } else { 21228c992abcSLaurent Vivier aSig &= LIT64(0xF800000000000000); 21238c992abcSLaurent Vivier aSig |= LIT64(0x0400000000000000); 21248c992abcSLaurent Vivier xsave = packFloatx80(aSign, aExp, aSig); /* F */ 21258c992abcSLaurent Vivier fp0 = a; 21268c992abcSLaurent Vivier fp1 = a; /* X */ 21278c992abcSLaurent Vivier fp2 = packFloatx80(0, one_exp, one_sig); /* 1 */ 21288c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, xsave, status); /* X*F */ 21298c992abcSLaurent Vivier fp0 = floatx80_sub(fp0, xsave, status); /* X-F */ 21308c992abcSLaurent Vivier fp1 = floatx80_add(fp1, fp2, status); /* 1 + X*F */ 21318c992abcSLaurent Vivier fp0 = floatx80_div(fp0, fp1, status); /* U = (X-F)/(1+X*F) */ 21328c992abcSLaurent Vivier 21338c992abcSLaurent Vivier tbl_index = compact; 21348c992abcSLaurent Vivier 21358c992abcSLaurent Vivier tbl_index &= 0x7FFF0000; 21368c992abcSLaurent Vivier tbl_index -= 0x3FFB0000; 21378c992abcSLaurent Vivier tbl_index >>= 1; 21388c992abcSLaurent Vivier tbl_index += compact & 0x00007800; 21398c992abcSLaurent Vivier tbl_index >>= 11; 21408c992abcSLaurent Vivier 21418c992abcSLaurent Vivier fp3 = atan_tbl[tbl_index]; 21428c992abcSLaurent Vivier 21438c992abcSLaurent Vivier fp3.high |= aSign ? 0x8000 : 0; /* ATAN(F) */ 21448c992abcSLaurent Vivier 21458c992abcSLaurent Vivier fp1 = floatx80_mul(fp0, fp0, status); /* V = U*U */ 21468c992abcSLaurent Vivier fp2 = float64_to_floatx80(make_float64(0xBFF6687E314987D8), 21478c992abcSLaurent Vivier status); /* A3 */ 21488c992abcSLaurent Vivier fp2 = floatx80_add(fp2, fp1, status); /* A3+V */ 21498c992abcSLaurent Vivier fp2 = floatx80_mul(fp2, fp1, status); /* V*(A3+V) */ 21508c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, fp0, status); /* U*V */ 21518c992abcSLaurent Vivier fp2 = floatx80_add(fp2, float64_to_floatx80( 21528c992abcSLaurent Vivier make_float64(0x4002AC6934A26DB3), status), 21538c992abcSLaurent Vivier status); /* A2+V*(A3+V) */ 21548c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, float64_to_floatx80( 21558c992abcSLaurent Vivier make_float64(0xBFC2476F4E1DA28E), status), 21568c992abcSLaurent Vivier status); /* A1+U*V */ 21578c992abcSLaurent Vivier fp1 = floatx80_mul(fp1, fp2, status); /* A1*U*V*(A2+V*(A3+V)) */ 21588c992abcSLaurent Vivier fp0 = floatx80_add(fp0, fp1, status); /* ATAN(U) */ 21598c992abcSLaurent Vivier 21608c992abcSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 21618c992abcSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 21628c992abcSLaurent Vivier 21638c992abcSLaurent Vivier a = floatx80_add(fp0, fp3, status); /* ATAN(X) */ 21648c992abcSLaurent Vivier 21658c992abcSLaurent Vivier float_raise(float_flag_inexact, status); 21668c992abcSLaurent Vivier 21678c992abcSLaurent Vivier return a; 21688c992abcSLaurent Vivier } 21698c992abcSLaurent Vivier } 2170*bc20b34eSLaurent Vivier 2171*bc20b34eSLaurent Vivier /*---------------------------------------------------------------------------- 2172*bc20b34eSLaurent Vivier | Arc sine 2173*bc20b34eSLaurent Vivier *----------------------------------------------------------------------------*/ 2174*bc20b34eSLaurent Vivier 2175*bc20b34eSLaurent Vivier floatx80 floatx80_asin(floatx80 a, float_status *status) 2176*bc20b34eSLaurent Vivier { 2177*bc20b34eSLaurent Vivier flag aSign; 2178*bc20b34eSLaurent Vivier int32_t aExp; 2179*bc20b34eSLaurent Vivier uint64_t aSig; 2180*bc20b34eSLaurent Vivier 2181*bc20b34eSLaurent Vivier int8_t user_rnd_mode, user_rnd_prec; 2182*bc20b34eSLaurent Vivier 2183*bc20b34eSLaurent Vivier int32_t compact; 2184*bc20b34eSLaurent Vivier floatx80 fp0, fp1, fp2, one; 2185*bc20b34eSLaurent Vivier 2186*bc20b34eSLaurent Vivier aSig = extractFloatx80Frac(a); 2187*bc20b34eSLaurent Vivier aExp = extractFloatx80Exp(a); 2188*bc20b34eSLaurent Vivier aSign = extractFloatx80Sign(a); 2189*bc20b34eSLaurent Vivier 2190*bc20b34eSLaurent Vivier if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) { 2191*bc20b34eSLaurent Vivier return propagateFloatx80NaNOneArg(a, status); 2192*bc20b34eSLaurent Vivier } 2193*bc20b34eSLaurent Vivier 2194*bc20b34eSLaurent Vivier if (aExp == 0 && aSig == 0) { 2195*bc20b34eSLaurent Vivier return packFloatx80(aSign, 0, 0); 2196*bc20b34eSLaurent Vivier } 2197*bc20b34eSLaurent Vivier 2198*bc20b34eSLaurent Vivier compact = floatx80_make_compact(aExp, aSig); 2199*bc20b34eSLaurent Vivier 2200*bc20b34eSLaurent Vivier if (compact >= 0x3FFF8000) { /* |X| >= 1 */ 2201*bc20b34eSLaurent Vivier if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */ 2202*bc20b34eSLaurent Vivier float_raise(float_flag_inexact, status); 2203*bc20b34eSLaurent Vivier a = packFloatx80(aSign, piby2_exp, pi_sig); 2204*bc20b34eSLaurent Vivier return floatx80_move(a, status); 2205*bc20b34eSLaurent Vivier } else { /* |X| > 1 */ 2206*bc20b34eSLaurent Vivier float_raise(float_flag_invalid, status); 2207*bc20b34eSLaurent Vivier return floatx80_default_nan(status); 2208*bc20b34eSLaurent Vivier } 2209*bc20b34eSLaurent Vivier 2210*bc20b34eSLaurent Vivier } /* |X| < 1 */ 2211*bc20b34eSLaurent Vivier 2212*bc20b34eSLaurent Vivier user_rnd_mode = status->float_rounding_mode; 2213*bc20b34eSLaurent Vivier user_rnd_prec = status->floatx80_rounding_precision; 2214*bc20b34eSLaurent Vivier status->float_rounding_mode = float_round_nearest_even; 2215*bc20b34eSLaurent Vivier status->floatx80_rounding_precision = 80; 2216*bc20b34eSLaurent Vivier 2217*bc20b34eSLaurent Vivier one = packFloatx80(0, one_exp, one_sig); 2218*bc20b34eSLaurent Vivier fp0 = a; 2219*bc20b34eSLaurent Vivier 2220*bc20b34eSLaurent Vivier fp1 = floatx80_sub(one, fp0, status); /* 1 - X */ 2221*bc20b34eSLaurent Vivier fp2 = floatx80_add(one, fp0, status); /* 1 + X */ 2222*bc20b34eSLaurent Vivier fp1 = floatx80_mul(fp2, fp1, status); /* (1+X)*(1-X) */ 2223*bc20b34eSLaurent Vivier fp1 = floatx80_sqrt(fp1, status); /* SQRT((1+X)*(1-X)) */ 2224*bc20b34eSLaurent Vivier fp0 = floatx80_div(fp0, fp1, status); /* X/SQRT((1+X)*(1-X)) */ 2225*bc20b34eSLaurent Vivier 2226*bc20b34eSLaurent Vivier status->float_rounding_mode = user_rnd_mode; 2227*bc20b34eSLaurent Vivier status->floatx80_rounding_precision = user_rnd_prec; 2228*bc20b34eSLaurent Vivier 2229*bc20b34eSLaurent Vivier a = floatx80_atan(fp0, status); /* ATAN(X/SQRT((1+X)*(1-X))) */ 2230*bc20b34eSLaurent Vivier 2231*bc20b34eSLaurent Vivier float_raise(float_flag_inexact, status); 2232*bc20b34eSLaurent Vivier 2233*bc20b34eSLaurent Vivier return a; 2234*bc20b34eSLaurent Vivier } 2235