libc/src/math/generic/atanf.cpp

   1 //===-- Single-precision atan 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/atanf.h"
  10 #include "inv_trigf_utils.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/FPUtil/rounding_mode.h"
  17 #include "src/__support/macros/config.h"
  18 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
  19
  20 namespace LIBC_NAMESPACE_DECL {
  21
  22 LLVM_LIBC_FUNCTION(float, atanf, (float x)) {
  23   using FPBits = typename fputil::FPBits<float>;
  24
  25   constexpr double FINAL_SIGN[2] = {1.0, -1.0};
  26   constexpr double SIGNED_PI_OVER_2[2] = {0x1.921fb54442d18p0,
  27                                           -0x1.921fb54442d18p0};
  28
  29   FPBits x_bits(x);
  30   Sign sign = x_bits.sign();
  31   x_bits.set_sign(Sign::POS);
  32   uint32_t x_abs = x_bits.uintval();
  33
  34   // x is inf or nan, |x| < 2^-4 or |x|= > 16.
  35   if (LIBC_UNLIKELY(x_abs <= 0x3d80'0000U || x_abs >= 0x4180'0000U)) {
  36     double x_d = static_cast<double>(x);
  37     double const_term = 0.0;
  38     if (LIBC_UNLIKELY(x_abs >= 0x4180'0000)) {
  39       // atan(+-Inf) = +-pi/2.
  40       if (x_bits.is_inf()) {
  41         volatile double sign_pi_over_2 = SIGNED_PI_OVER_2[sign.is_neg()];
  42         return static_cast<float>(sign_pi_over_2);
  43       }
  44       if (x_bits.is_nan())
  45         return x;
  46       // x >= 16
  47       x_d = -1.0 / x_d;
  48       const_term = SIGNED_PI_OVER_2[sign.is_neg()];
  49     }
  50     // 0 <= x < 1/16;
  51     if (LIBC_UNLIKELY(x_bits.is_zero()))
  52       return x;
  53     // x <= 2^-12;
  54     if (LIBC_UNLIKELY(x_abs < 0x3980'0000)) {
  55 #if defined(LIBC_TARGET_CPU_HAS_FMA)
  56       return fputil::multiply_add(x, -0x1.0p-25f, x);
  57 #else
  58       double x_d = static_cast<double>(x);
  59       return static_cast<float>(fputil::multiply_add(x_d, -0x1.0p-25, x_d));
  60 #endif // LIBC_TARGET_CPU_HAS_FMA
  61     }
  62     // Use Taylor polynomial:
  63     //   atan(x) ~ x * (1 - x^2 / 3 + x^4 / 5 - x^6 / 7 + x^8 / 9 - x^10 / 11).
  64     constexpr double ATAN_TAYLOR[6] = {
  65         0x1.0000000000000p+0,  -0x1.5555555555555p-2, 0x1.999999999999ap-3,
  66         -0x1.2492492492492p-3, 0x1.c71c71c71c71cp-4,  -0x1.745d1745d1746p-4,
  67     };
  68     double x2 = x_d * x_d;
  69     double x4 = x2 * x2;
  70     double c0 = fputil::multiply_add(x2, ATAN_TAYLOR[1], ATAN_TAYLOR[0]);
  71     double c1 = fputil::multiply_add(x2, ATAN_TAYLOR[3], ATAN_TAYLOR[2]);
  72     double c2 = fputil::multiply_add(x2, ATAN_TAYLOR[5], ATAN_TAYLOR[4]);
  73     double p = fputil::polyeval(x4, c0, c1, c2);
  74     double r = fputil::multiply_add(x_d, p, const_term);
  75     return static_cast<float>(r);
  76   }
  77
  78   // Range reduction steps:
  79   // 1)  atan(x) = sign(x) * atan(|x|)
  80   // 2)  If |x| > 1, atan(|x|) = pi/2 - atan(1/|x|)
  81   // 3)  For 1/16 < x <= 1, we find k such that: |x - k/16| <= 1/32.
  82   // 4)  Then we use polynomial approximation:
  83   //   atan(x) ~ atan((k/16) + (x - (k/16)) * Q(x - k/16)
  84   //           = P(x - k/16)
  85   double x_d, const_term, final_sign;
  86   int idx;
  87
  88   if (x_abs > 0x3f80'0000U) {
  89     // |x| > 1, we need to invert x, so we will perform range reduction in
  90     // double precision.
  91     x_d = 1.0 / static_cast<double>(x_bits.get_val());
  92     double k_d = fputil::nearest_integer(x_d * 0x1.0p4);
  93     x_d = fputil::multiply_add(k_d, -0x1.0p-4, x_d);
  94     idx = static_cast<int>(k_d);
  95     final_sign = FINAL_SIGN[sign.is_pos()];
  96     // Adjust constant term of the polynomial by +- pi/2.
  97     const_term = fputil::multiply_add(final_sign, ATAN_COEFFS[idx][0],
  98                                       SIGNED_PI_OVER_2[sign.is_neg()]);
  99   } else {
 100     // Exceptional value:
 101     if (LIBC_UNLIKELY(x_abs == 0x3d8d'6b23U)) { // |x| = 0x1.1ad646p-4
 102       return sign.is_pos() ? fputil::round_result_slightly_down(0x1.1a6386p-4f)
 103                            : fputil::round_result_slightly_up(-0x1.1a6386p-4f);
 104     }
 105     // Perform range reduction in single precision.
 106     float x_f = x_bits.get_val();
 107     float k_f = fputil::nearest_integer(x_f * 0x1.0p4f);
 108     x_f = fputil::multiply_add(k_f, -0x1.0p-4f, x_f);
 109     x_d = static_cast<double>(x_f);
 110     idx = static_cast<int>(k_f);
 111     final_sign = FINAL_SIGN[sign.is_neg()];
 112     const_term = final_sign * ATAN_COEFFS[idx][0];
 113   }
 114
 115   double p = atan_eval(x_d, idx);
 116   double r = fputil::multiply_add(final_sign * x_d, p, const_term);
 117
 118   return static_cast<float>(r);
 119 }
 120
 121 } // namespace LIBC_NAMESPACE_DECL