1 //===-- Single-precision tan 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/tanf.h" 10 #include "sincosf_utils.h" 11 #include "src/__support/FPUtil/FEnvImpl.h" 12 #include "src/__support/FPUtil/FPBits.h" 13 #include "src/__support/FPUtil/PolyEval.h" 14 #include "src/__support/FPUtil/except_value_utils.h" 15 #include "src/__support/FPUtil/multiply_add.h" 16 #include "src/__support/FPUtil/nearest_integer.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 22 namespace LIBC_NAMESPACE_DECL { 23 24 // Exceptional cases for tanf. 25 constexpr size_t N_EXCEPTS = 6; 26 27 constexpr fputil::ExceptValues<float, N_EXCEPTS> TANF_EXCEPTS{{ 28 // (inputs, RZ output, RU offset, RD offset, RN offset) 29 // x = 0x1.ada6aap27, tan(x) = 0x1.e80304p-3 (RZ) 30 {0x4d56d355, 0x3e740182, 1, 0, 0}, 31 // x = 0x1.862064p33, tan(x) = -0x1.8dee56p-3 (RZ) 32 {0x50431032, 0xbe46f72b, 0, 1, 1}, 33 // x = 0x1.af61dap48, tan(x) = 0x1.60d1c6p-2 (RZ) 34 {0x57d7b0ed, 0x3eb068e3, 1, 0, 1}, 35 // x = 0x1.0088bcp52, tan(x) = 0x1.ca1edp0 (RZ) 36 {0x5980445e, 0x3fe50f68, 1, 0, 0}, 37 // x = 0x1.f90dfcp72, tan(x) = 0x1.597f9cp-1 (RZ) 38 {0x63fc86fe, 0x3f2cbfce, 1, 0, 0}, 39 // x = 0x1.a6ce12p86, tan(x) = -0x1.c5612ep-1 (RZ) 40 {0x6ad36709, 0xbf62b097, 0, 1, 0}, 41 }}; 42 43 LLVM_LIBC_FUNCTION(float, tanf, (float x)) { 44 using FPBits = typename fputil::FPBits<float>; 45 FPBits xbits(x); 46 bool x_sign = xbits.uintval() >> 31; 47 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU; 48 49 // |x| < pi/32 50 if (LIBC_UNLIKELY(x_abs <= 0x3dc9'0fdbU)) { 51 double xd = static_cast<double>(x); 52 53 // |x| < 0x1.0p-12f 54 if (LIBC_UNLIKELY(x_abs < 0x3980'0000U)) { 55 if (LIBC_UNLIKELY(x_abs == 0U)) { 56 // For signed zeros. 57 return x; 58 } 59 // When |x| < 2^-12, the relative error of the approximation tan(x) ~ x 60 // is: 61 // |tan(x) - x| / |tan(x)| < |x^3| / (3|x|) 62 // = x^2 / 3 63 // < 2^-25 64 // < epsilon(1)/2. 65 // So the correctly rounded values of tan(x) are: 66 // = x + sign(x)*eps(x) if rounding mode = FE_UPWARD and x is positive, 67 // or (rounding mode = FE_DOWNWARD and x is 68 // negative), 69 // = x otherwise. 70 // To simplify the rounding decision and make it more efficient, we use 71 // fma(x, 2^-25, x) instead. 72 // Note: to use the formula x + 2^-25*x to decide the correct rounding, we 73 // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when 74 // |x| < 2^-125. For targets without FMA instructions, we simply use 75 // double for intermediate results as it is more efficient than using an 76 // emulated version of FMA. 77 #if defined(LIBC_TARGET_CPU_HAS_FMA) 78 return fputil::multiply_add(x, 0x1.0p-25f, x); 79 #else 80 return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd)); 81 #endif // LIBC_TARGET_CPU_HAS_FMA 82 } 83 84 // |x| < pi/32 85 double xsq = xd * xd; 86 87 // Degree-9 minimax odd polynomial of tan(x) generated by Sollya with: 88 // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6, 8|], [|1, D...|], [0, pi/32]); 89 double result = 90 fputil::polyeval(xsq, 1.0, 0x1.555555553d022p-2, 0x1.111111ce442c1p-3, 91 0x1.ba180a6bbdecdp-5, 0x1.69c0a88a0b71fp-6); 92 return static_cast<float>(xd * result); 93 } 94 95 // Check for exceptional values 96 if (LIBC_UNLIKELY(x_abs == 0x3f8a1f62U)) { 97 // |x| = 0x1.143ec4p0 98 float sign = x_sign ? -1.0f : 1.0f; 99 100 // volatile is used to prevent compiler (gcc) from optimizing the 101 // computation, making the results incorrect in different rounding modes. 102 volatile float tmp = 0x1.ddf9f4p0f; 103 tmp = fputil::multiply_add(sign, tmp, sign * 0x1.1p-24f); 104 105 return tmp; 106 } 107 108 // |x| > 0x1.ada6a8p+27f 109 if (LIBC_UNLIKELY(x_abs > 0x4d56'd354U)) { 110 // Inf or NaN 111 if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) { 112 if (x_abs == 0x7f80'0000U) { 113 fputil::set_errno_if_required(EDOM); 114 fputil::raise_except_if_required(FE_INVALID); 115 } 116 return x + FPBits::quiet_nan().get_val(); 117 } 118 // Other large exceptional values 119 if (auto r = TANF_EXCEPTS.lookup_odd(x_abs, x_sign); 120 LIBC_UNLIKELY(r.has_value())) 121 return r.value(); 122 } 123 124 // For |x| >= pi/32, we use the definition of tan(x) function: 125 // tan(x) = sin(x) / cos(x) 126 // The we follow the same computations of sin(x) and cos(x) as sinf, cosf, 127 // and sincosf. 128 129 double xd = static_cast<double>(x); 130 double sin_k, cos_k, sin_y, cosm1_y; 131 132 sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y); 133 // tan(x) = sin(x) / cos(x) 134 // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k) 135 using fputil::multiply_add; 136 return static_cast<float>( 137 multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) / 138 multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k))); 139 } 140 141 } // namespace LIBC_NAMESPACE_DECL 142