153 lines · cpp
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-exception6//7//===----------------------------------------------------------------------===//8 9#include "src/math/tanf.h"10#include "src/__support/FPUtil/FEnvImpl.h"11#include "src/__support/FPUtil/FPBits.h"12#include "src/__support/FPUtil/PolyEval.h"13#include "src/__support/FPUtil/except_value_utils.h"14#include "src/__support/FPUtil/multiply_add.h"15#include "src/__support/FPUtil/nearest_integer.h"16#include "src/__support/common.h"17#include "src/__support/macros/config.h"18#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY19#include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA20#include "src/__support/math/sincosf_utils.h"21 22namespace LIBC_NAMESPACE_DECL {23 24#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS25// Exceptional cases for tanf.26constexpr size_t N_EXCEPTS = 6;27 28constexpr fputil::ExceptValues<float, N_EXCEPTS> TANF_EXCEPTS{{29 // (inputs, RZ output, RU offset, RD offset, RN offset)30 // x = 0x1.ada6aap27, tan(x) = 0x1.e80304p-3 (RZ)31 {0x4d56d355, 0x3e740182, 1, 0, 0},32 // x = 0x1.862064p33, tan(x) = -0x1.8dee56p-3 (RZ)33 {0x50431032, 0xbe46f72b, 0, 1, 1},34 // x = 0x1.af61dap48, tan(x) = 0x1.60d1c6p-2 (RZ)35 {0x57d7b0ed, 0x3eb068e3, 1, 0, 1},36 // x = 0x1.0088bcp52, tan(x) = 0x1.ca1edp0 (RZ)37 {0x5980445e, 0x3fe50f68, 1, 0, 0},38 // x = 0x1.f90dfcp72, tan(x) = 0x1.597f9cp-1 (RZ)39 {0x63fc86fe, 0x3f2cbfce, 1, 0, 0},40 // x = 0x1.a6ce12p86, tan(x) = -0x1.c5612ep-1 (RZ)41 {0x6ad36709, 0xbf62b097, 0, 1, 0},42}};43#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS44 45LLVM_LIBC_FUNCTION(float, tanf, (float x)) {46 using FPBits = typename fputil::FPBits<float>;47 FPBits xbits(x);48 uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;49 50 // |x| < pi/3251 if (LIBC_UNLIKELY(x_abs <= 0x3dc9'0fdbU)) {52 double xd = static_cast<double>(x);53 54 // |x| < 0x1.0p-12f55 if (LIBC_UNLIKELY(x_abs < 0x3980'0000U)) {56 if (LIBC_UNLIKELY(x_abs == 0U)) {57 // For signed zeros.58 return x;59 }60 // When |x| < 2^-12, the relative error of the approximation tan(x) ~ x61 // is:62 // |tan(x) - x| / |tan(x)| < |x^3| / (3|x|)63 // = x^2 / 364 // < 2^-2565 // < epsilon(1)/2.66 // So the correctly rounded values of tan(x) are:67 // = x + sign(x)*eps(x) if rounding mode = FE_UPWARD and x is positive,68 // or (rounding mode = FE_DOWNWARD and x is69 // negative),70 // = x otherwise.71 // To simplify the rounding decision and make it more efficient, we use72 // fma(x, 2^-25, x) instead.73 // Note: to use the formula x + 2^-25*x to decide the correct rounding, we74 // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when75 // |x| < 2^-125. For targets without FMA instructions, we simply use76 // double for intermediate results as it is more efficient than using an77 // emulated version of FMA.78#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)79 return fputil::multiply_add(x, 0x1.0p-25f, x);80#else81 return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));82#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT83 }84 85 // |x| < pi/3286 double xsq = xd * xd;87 88 // Degree-9 minimax odd polynomial of tan(x) generated by Sollya with:89 // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6, 8|], [|1, D...|], [0, pi/32]);90 double result =91 fputil::polyeval(xsq, 1.0, 0x1.555555553d022p-2, 0x1.111111ce442c1p-3,92 0x1.ba180a6bbdecdp-5, 0x1.69c0a88a0b71fp-6);93 return static_cast<float>(xd * result);94 }95 96#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS97 bool x_sign = xbits.uintval() >> 31;98 // Check for exceptional values99 if (LIBC_UNLIKELY(x_abs == 0x3f8a1f62U)) {100 // |x| = 0x1.143ec4p0101 float sign = x_sign ? -1.0f : 1.0f;102 103 // volatile is used to prevent compiler (gcc) from optimizing the104 // computation, making the results incorrect in different rounding modes.105 volatile float tmp = 0x1.ddf9f4p0f;106 tmp = fputil::multiply_add(sign, tmp, sign * 0x1.1p-24f);107 108 return tmp;109 }110#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS111 112 // |x| > 0x1.ada6a8p+27f113 if (LIBC_UNLIKELY(x_abs > 0x4d56'd354U)) {114 // Inf or NaN115 if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {116 if (xbits.is_signaling_nan()) {117 fputil::raise_except_if_required(FE_INVALID);118 return FPBits::quiet_nan().get_val();119 }120 121 if (x_abs == 0x7f80'0000U) {122 fputil::set_errno_if_required(EDOM);123 fputil::raise_except_if_required(FE_INVALID);124 }125 return x + FPBits::quiet_nan().get_val();126 }127#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS128 // Other large exceptional values129 if (auto r = TANF_EXCEPTS.lookup_odd(x_abs, x_sign);130 LIBC_UNLIKELY(r.has_value()))131 return r.value();132#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS133 }134 135 // For |x| >= pi/32, we use the definition of tan(x) function:136 // tan(x) = sin(x) / cos(x)137 // The we follow the same computations of sin(x) and cos(x) as sinf, cosf,138 // and sincosf.139 140 double xd = static_cast<double>(x);141 double sin_k, cos_k, sin_y, cosm1_y;142 143 sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);144 // tan(x) = sin(x) / cos(x)145 // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)146 using fputil::multiply_add;147 return static_cast<float>(148 multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) /149 multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k)));150}151 152} // namespace LIBC_NAMESPACE_DECL153