1 //===-- Double-precision cos 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/cos.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/except_value_utils.h" 16 #include "src/__support/common.h" 17 #include "src/__support/macros/config.h" 18 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY 19 #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA 20 #include "src/math/generic/sincos_eval.h" 21 22 #ifdef LIBC_TARGET_CPU_HAS_FMA 23 #include "range_reduction_double_fma.h" 24 25 using LIBC_NAMESPACE::fma::FAST_PASS_EXPONENT; 26 using LIBC_NAMESPACE::fma::ONE_TWENTY_EIGHT_OVER_PI; 27 using LIBC_NAMESPACE::fma::range_reduction_small; 28 using LIBC_NAMESPACE::fma::SIN_K_PI_OVER_128; 29 30 LIBC_INLINE constexpr bool NO_FMA = false; 31 #else 32 #include "range_reduction_double_nofma.h" 33 34 using LIBC_NAMESPACE::nofma::FAST_PASS_EXPONENT; 35 using LIBC_NAMESPACE::nofma::ONE_TWENTY_EIGHT_OVER_PI; 36 using LIBC_NAMESPACE::nofma::range_reduction_small; 37 using LIBC_NAMESPACE::nofma::SIN_K_PI_OVER_128; 38 39 LIBC_INLINE constexpr bool NO_FMA = true; 40 #endif // LIBC_TARGET_CPU_HAS_FMA 41 42 // TODO: We might be able to improve the performance of large range reduction of 43 // non-FMA targets further by operating directly on 25-bit chunks of 128/pi and 44 // pre-split SIN_K_PI_OVER_128, but that might double the memory footprint of 45 // those lookup table. 46 #include "range_reduction_double_common.h" 47 48 #if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0) 49 #define LIBC_MATH_COS_SKIP_ACCURATE_PASS 50 #endif 51 52 namespace LIBC_NAMESPACE_DECL { 53 54 using DoubleDouble = fputil::DoubleDouble; 55 using Float128 = typename fputil::DyadicFloat<128>; 56 57 LLVM_LIBC_FUNCTION(double, cos, (double x)) { 58 using FPBits = typename fputil::FPBits<double>; 59 FPBits xbits(x); 60 61 uint16_t x_e = xbits.get_biased_exponent(); 62 63 DoubleDouble y; 64 unsigned k; 65 generic::LargeRangeReduction<NO_FMA> range_reduction_large{}; 66 67 // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA) 68 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) { 69 // |x| < 2^-27 70 if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) { 71 // Signed zeros. 72 if (LIBC_UNLIKELY(x == 0.0)) 73 return 1.0; 74 75 // For |x| < 2^-27, |cos(x) - 1| < |x|^2/2 < 2^-54 = ulp(1 - 2^-53)/2. 76 return fputil::round_result_slightly_down(1.0); 77 } 78 79 // // Small range reduction. 80 k = range_reduction_small(x, y); 81 } else { 82 // Inf or NaN 83 if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) { 84 // sin(+-Inf) = NaN 85 if (xbits.get_mantissa() == 0) { 86 fputil::set_errno_if_required(EDOM); 87 fputil::raise_except_if_required(FE_INVALID); 88 } 89 return x + FPBits::quiet_nan().get_val(); 90 } 91 92 // Large range reduction. 93 k = range_reduction_large.compute_high_part(x); 94 y = range_reduction_large.fast(); 95 } 96 97 DoubleDouble sin_y, cos_y; 98 99 generic::sincos_eval(y, sin_y, cos_y); 100 101 // Look up sin(k * pi/128) and cos(k * pi/128) 102 // Memory saving versions: 103 104 // Use 128-entry table instead: 105 // DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 127]; 106 // uint64_t sin_s = static_cast<uint64_t>((k + 128) & 128) << (63 - 7); 107 // sin_k.hi = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); 108 // sin_k.lo = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val(); 109 // DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 127]; 110 // uint64_t cos_s = static_cast<uint64_t>((k + 64) & 128) << (63 - 7); 111 // cos_k.hi = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); 112 // cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val(); 113 114 // Use 64-entry table instead: 115 // auto get_idx_dd = [](unsigned kk) -> DoubleDouble { 116 // unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); 117 // DoubleDouble ans = SIN_K_PI_OVER_128[idx]; 118 // if (kk & 128) { 119 // ans.hi = -ans.hi; 120 // ans.lo = -ans.lo; 121 // } 122 // return ans; 123 // }; 124 // DoubleDouble sin_k = get_idx_dd(k + 128); 125 // DoubleDouble cos_k = get_idx_dd(k + 64); 126 127 // Fast look up version, but needs 256-entry table. 128 // -sin(k * pi/128) = sin((k + 128) * pi/128) 129 // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). 130 DoubleDouble msin_k = SIN_K_PI_OVER_128[(k + 128) & 255]; 131 DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255]; 132 133 // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128). 134 // So k is an integer and -pi / 256 <= y <= pi / 256. 135 // Then cos(x) = cos((k * pi/128 + y) 136 // = cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128) 137 DoubleDouble cos_k_cos_y = fputil::quick_mult<NO_FMA>(cos_y, cos_k); 138 DoubleDouble msin_k_sin_y = fputil::quick_mult<NO_FMA>(sin_y, msin_k); 139 140 DoubleDouble rr = fputil::exact_add<false>(cos_k_cos_y.hi, msin_k_sin_y.hi); 141 rr.lo += msin_k_sin_y.lo + cos_k_cos_y.lo; 142 143 #ifdef LIBC_MATH_COS_SKIP_ACCURATE_PASS 144 return rr.hi + rr.lo; 145 #else 146 147 // Accurate test and pass for correctly rounded implementation. 148 #ifdef LIBC_TARGET_CPU_HAS_FMA 149 constexpr double ERR = 0x1.0p-70; 150 #else 151 // TODO: Improve non-FMA fast pass accuracy. 152 constexpr double ERR = 0x1.0p-66; 153 #endif // LIBC_TARGET_CPU_HAS_FMA 154 155 double rlp = rr.lo + ERR; 156 double rlm = rr.lo - ERR; 157 158 double r_upper = rr.hi + rlp; // (rr.lo + ERR); 159 double r_lower = rr.hi + rlm; // (rr.lo - ERR); 160 161 // Ziv's rounding test. 162 if (LIBC_LIKELY(r_upper == r_lower)) 163 return r_upper; 164 165 Float128 u_f128, sin_u, cos_u; 166 if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) 167 u_f128 = generic::range_reduction_small_f128(x); 168 else 169 u_f128 = range_reduction_large.accurate(); 170 171 generic::sincos_eval(u_f128, sin_u, cos_u); 172 173 auto get_sin_k = [](unsigned kk) -> Float128 { 174 unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63); 175 Float128 ans = generic::SIN_K_PI_OVER_128_F128[idx]; 176 if (kk & 128) 177 ans.sign = Sign::NEG; 178 return ans; 179 }; 180 181 // -sin(k * pi/128) = sin((k + 128) * pi/128) 182 // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128). 183 Float128 msin_k_f128 = get_sin_k(k + 128); 184 Float128 cos_k_f128 = get_sin_k(k + 64); 185 186 // cos(x) = cos((k * pi/128 + u) 187 // = cos(u) * cos(k*pi/128) - sin(u) * sin(k*pi/128) 188 Float128 r = fputil::quick_add(fputil::quick_mul(cos_k_f128, cos_u), 189 fputil::quick_mul(msin_k_f128, sin_u)); 190 191 // TODO: Add assertion if Ziv's accuracy tests fail in debug mode. 192 // https://github.com/llvm/llvm-project/issues/96452. 193 194 return static_cast<double>(r); 195 #endif // !LIBC_MATH_COS_SKIP_ACCURATE_PASS 196 } 197 198 } // namespace LIBC_NAMESPACE_DECL 199