libc/src/math/generic/tanf.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-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/optimization.h"            // LIBC_UNLIKELY
  19 #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
  20
  21 #include <errno.h>
  22
  23 namespace __llvm_libc {
  24
  25 // Exceptional cases for tanf.
  26 constexpr size_t N_EXCEPTS = 6;
  27
  28 constexpr 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
  44 LLVM_LIBC_FUNCTION(float, tanf, (float x)) {
  45   using FPBits = typename fputil::FPBits<float>;
  46   FPBits xbits(x);
  47   bool x_sign = xbits.uintval() >> 31;
  48   uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
  49
  50   // |x| < pi/32
  51   if (LIBC_UNLIKELY(x_abs <= 0x3dc9'0fdbU)) {
  52     double xd = static_cast<double>(x);
  53
  54     // |x| < 0x1.0p-12f
  55     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) ~ x
  61       // is:
  62       //   |tan(x) - x| / |tan(x)| < |x^3| / (3|x|)
  63       //                           = x^2 / 3
  64       //                           < 2^-25
  65       //                           < 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 is
  69       //                        negative),
  70       //   = x otherwise.
  71       // To simplify the rounding decision and make it more efficient, we use
  72       //   fma(x, 2^-25, x) instead.
  73       // Note: to use the formula x + 2^-25*x to decide the correct rounding, we
  74       // do need fma(x, 2^-25, x) to prevent underflow caused by 2^-25*x when
  75       // |x| < 2^-125. For targets without FMA instructions, we simply use
  76       // double for intermediate results as it is more efficient than using an
  77       // emulated version of FMA.
  78 #if defined(LIBC_TARGET_CPU_HAS_FMA)
  79       return fputil::multiply_add(x, 0x1.0p-25f, x);
  80 #else
  81       return static_cast<float>(fputil::multiply_add(xd, 0x1.0p-25, xd));
  82 #endif // LIBC_TARGET_CPU_HAS_FMA
  83     }
  84
  85     // |x| < pi/32
  86     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   // Check for exceptional values
  97   if (LIBC_UNLIKELY(x_abs == 0x3f8a1f62U)) {
  98     // |x| = 0x1.143ec4p0
  99     float sign = x_sign ? -1.0f : 1.0f;
 100
 101     // volatile is used to prevent compiler (gcc) from optimizing the
 102     // computation, making the results incorrect in different rounding modes.
 103     volatile float tmp = 0x1.ddf9f4p0f;
 104     tmp = fputil::multiply_add(sign, tmp, sign * 0x1.1p-24f);
 105
 106     return tmp;
 107   }
 108
 109   // |x| > 0x1.ada6a8p+27f
 110   if (LIBC_UNLIKELY(x_abs > 0x4d56'd354U)) {
 111     // Inf or NaN
 112     if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
 113       if (x_abs == 0x7f80'0000U) {
 114         fputil::set_errno_if_required(EDOM);
 115         fputil::raise_except_if_required(FE_INVALID);
 116       }
 117       return x + FPBits::build_quiet_nan(0);
 118     }
 119     // Other large exceptional values
 120     if (auto r = TANF_EXCEPTS.lookup_odd(x_abs, x_sign);
 121         LIBC_UNLIKELY(r.has_value()))
 122       return r.value();
 123   }
 124
 125   // For |x| >= pi/32, we use the definition of tan(x) function:
 126   //   tan(x) = sin(x) / cos(x)
 127   // The we follow the same computations of sin(x) and cos(x) as sinf, cosf,
 128   // and sincosf.
 129
 130   double xd = static_cast<double>(x);
 131   double sin_k, cos_k, sin_y, cosm1_y;
 132
 133   sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
 134   // tan(x) = sin(x) / cos(x)
 135   //        = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
 136   using fputil::multiply_add;
 137   return static_cast<float>(
 138       multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) /
 139       multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k)));
 140 }
 141
 142 } // namespace __llvm_libc