1 //===-- Double-precision sin function -------------------------------------===// 2 // 3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. 4 // See https://llvm.org/LICENSE.txt for license information. 5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception 6 // 7 //===----------------------------------------------------------------------===// 8 9 #include "src/math/sin.h" 10 #include "hdr/errno_macros.h" 11 #include "src/__support/FPUtil/FEnvImpl.h" 12 #include "src/__support/FPUtil/FPBits.h" 13 #include "src/__support/FPUtil/double_double.h" 14 #include "src/__support/FPUtil/dyadic_float.h" 15 #include "src/__support/FPUtil/multiply_add.h" 16 #include "src/__support/FPUtil/rounding_mode.h" 17 #include "src/__support/common.h" 18 #include "src/__support/macros/config.h" 19 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY 20 #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA 21 #include "src/math/generic/range_reduction_double_common.h" 22 #include "src/math/generic/sincos_eval.h" 23 24 #ifdef LIBC_TARGET_CPU_HAS_FMA 25 #include "range_reduction_double_fma.h" 26 #else 27 #include "range_reduction_double_nofma.h" 28 #endif // LIBC_TARGET_CPU_HAS_FMA 29 30 namespace LIBC_NAMESPACE_DECL { 31 32 using DoubleDouble = fputil::DoubleDouble; 33 using Float128 = typename fputil::DyadicFloat<128>; 34 35 LLVM_LIBC_FUNCTION(double, sin, (double x)) { 36 using FPBits = typename fputil::FPBits<double>; 37 FPBits xbits(x); 38 39 uint16_t x_e = xbits.get_biased_exponent(); 40 41 DoubleDouble y; 42 unsigned k; 43 LargeRangeReduction range_reduction_large{}; 44 45 // |x| < 2^16 46 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { 47 // |x| < 2^-7 48 if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) { 49 // |x| < 2^-26, |sin(x) - x| < ulp(x)/2. 50 if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 26)) { 51 // Signed zeros. 52 if (LIBC_UNLIKELY(x == 0.0)) 53 return x + x; // Make sure it works with FTZ/DAZ. 54 55 #ifdef LIBC_TARGET_CPU_HAS_FMA 56 return fputil::multiply_add(x, -0x1.0p-54, x); 57 #else 58 if (LIBC_UNLIKELY(x_e < 4)) { 59 int rounding_mode = fputil::quick_get_round(); 60 if (rounding_mode == FE_TOWARDZERO || 61 (xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) || 62 (xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD)) 63 return FPBits(xbits.uintval() - 1).get_val(); 64 } 65 return fputil::multiply_add(x, -0x1.0p-54, x); 66 #endif // LIBC_TARGET_CPU_HAS_FMA 67 } 68 // No range reduction needed. 69 k = 0; 70 y.lo = 0.0; 71 y.hi = x; 72 } else { 73 // Small range reduction. 74 k = range_reduction_small(x, y); 75 } 76 } else { 77 // Inf or NaN 78 if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) { 79 // sin(+-Inf) = NaN 80 if (xbits.get_mantissa() == 0) { 81 fputil::set_errno_if_required(EDOM); 82 fputil::raise_except_if_required(FE_INVALID); 83 } 84 return x + FPBits::quiet_nan().get_val(); 85 } 86 87 // Large range reduction. 88 k = range_reduction_large.fast(x, y); 89 } 90 91 DoubleDouble sin_y, cos_y; 92 93 [[maybe_unused]] double err = generic::sincos_eval(y, sin_y, cos_y); 94 95 // Look up sin(k * pi/128) and cos(k * pi/128) 96 #ifdef LIBC_MATH_HAS_SMALL_TABLES 97 // Memory saving versions. Use 65-entry table. __anonf1e1dc4f0102(unsigned kk) 98 auto get_idx_dd = [](unsigned kk) -> DoubleDouble { 99 unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); 100 DoubleDouble ans = SIN_K_PI_OVER_128[idx]; 101 if (kk & 128) { 102 ans.hi = -ans.hi; 103 ans.lo = -ans.lo; 104 } 105 return ans; 106 }; 107 DoubleDouble sin_k = get_idx_dd(k); 108 DoubleDouble cos_k = get_idx_dd(k + 64); 109 #else 110 // Fast look up version, but needs 256-entry table. 111 // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). 112 DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 255]; 113 DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255]; 114 #endif 115 116 // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128). 117 // So k is an integer and -pi / 256 <= y <= pi / 256. 118 // Then sin(x) = sin((k * pi/128 + y) 119 // = sin(y) * cos(k*pi/128) + cos(y) * sin(k*pi/128) 120 DoubleDouble sin_k_cos_y = fputil::quick_mult(cos_y, sin_k); 121 DoubleDouble cos_k_sin_y = fputil::quick_mult(sin_y, cos_k); 122 123 DoubleDouble rr = fputil::exact_add<false>(sin_k_cos_y.hi, cos_k_sin_y.hi); 124 rr.lo += sin_k_cos_y.lo + cos_k_sin_y.lo; 125 126 #ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS 127 return rr.hi + rr.lo; 128 #else 129 // Accurate test and pass for correctly rounded implementation. 130 131 double rlp = rr.lo + err; 132 double rlm = rr.lo - err; 133 134 double r_upper = rr.hi + rlp; // (rr.lo + ERR); 135 double r_lower = rr.hi + rlm; // (rr.lo - ERR); 136 137 // Ziv's rounding test. 138 if (LIBC_LIKELY(r_upper == r_lower)) 139 return r_upper; 140 141 Float128 u_f128, sin_u, cos_u; 142 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) 143 u_f128 = range_reduction_small_f128(x); 144 else 145 u_f128 = range_reduction_large.accurate(); 146 147 generic::sincos_eval(u_f128, sin_u, cos_u); 148 __anonf1e1dc4f0202(unsigned kk) 149 auto get_sin_k = [](unsigned kk) -> Float128 { 150 unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); 151 Float128 ans = SIN_K_PI_OVER_128_F128[idx]; 152 if (kk & 128) 153 ans.sign = Sign::NEG; 154 return ans; 155 }; 156 157 // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). 158 Float128 sin_k_f128 = get_sin_k(k); 159 Float128 cos_k_f128 = get_sin_k(k + 64); 160 161 // sin(x) = sin((k * pi/128 + u) 162 // = sin(u) * cos(k*pi/128) + cos(u) * sin(k*pi/128) 163 Float128 r = fputil::quick_add(fputil::quick_mul(sin_k_f128, cos_u), 164 fputil::quick_mul(cos_k_f128, sin_u)); 165 166 // TODO: Add assertion if Ziv's accuracy tests fail in debug mode. 167 // https://github.com/llvm/llvm-project/issues/96452. 168 169 return static_cast<double>(r); 170 #endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS 171 } 172 173 } // namespace LIBC_NAMESPACE_DECL 174