libc/src/math/generic/inv_trigf_utils.h

   1 //===-- Single-precision general inverse trigonometric functions ----------===//
   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 #ifndef LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H
  10 #define LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H
  11
  12 #include "math_utils.h"
  13 #include "src/__support/FPUtil/FEnvImpl.h"
  14 #include "src/__support/FPUtil/FPBits.h"
  15 #include "src/__support/FPUtil/PolyEval.h"
  16 #include "src/__support/FPUtil/nearest_integer.h"
  17 #include "src/__support/common.h"
  18
  19 #include <errno.h>
  20
  21 namespace __llvm_libc {
  22
  23 // PI and PI / 2
  24 constexpr double M_MATH_PI = 0x1.921fb54442d18p+1;
  25 constexpr double M_MATH_PI_2 = 0x1.921fb54442d18p+0;
  26
  27 // atan table size
  28 constexpr int ATAN_T_BITS = 4;
  29 constexpr int ATAN_T_SIZE = 1 << ATAN_T_BITS;
  30
  31 // N[Table[ArcTan[x], {x, 1/8, 8/8, 1/8}], 40]
  32 extern const double ATAN_T[ATAN_T_SIZE];
  33 extern const double ATAN_K[5];
  34
  35 // The main idea of the function is to use formula
  36 // atan(u) + atan(v) = atan((u+v)/(1-uv))
  37
  38 // x should be positive, normal finite value
  39 LIBC_INLINE double atan_eval(double x) {
  40   using FPB = fputil::FPBits<double>;
  41   // Added some small value to umin and umax mantissa to avoid possible rounding
  42   // errors.
  43   FPB::UIntType umin =
  44       FPB::create_value(false, FPB::EXPONENT_BIAS - ATAN_T_BITS - 1,
  45                         0x100000000000UL)
  46           .uintval();
  47   FPB::UIntType umax =
  48       FPB::create_value(false, FPB::EXPONENT_BIAS + ATAN_T_BITS,
  49                         0xF000000000000UL)
  50           .uintval();
  51
  52   FPB bs(x);
  53   bool sign = bs.get_sign();
  54   auto x_abs = bs.uintval() & FPB::FloatProp::EXP_MANT_MASK;
  55
  56   if (x_abs <= umin) {
  57     double pe = __llvm_libc::fputil::polyeval(x * x, 0.0, ATAN_K[1], ATAN_K[2],
  58                                               ATAN_K[3], ATAN_K[4]);
  59     return fputil::multiply_add(pe, x, x);
  60   }
  61
  62   if (x_abs >= umax) {
  63     double one_over_x_m = -1.0 / x;
  64     double one_over_x2 = one_over_x_m * one_over_x_m;
  65     double pe = __llvm_libc::fputil::polyeval(one_over_x2, ATAN_K[0], ATAN_K[1],
  66                                               ATAN_K[2], ATAN_K[3]);
  67     return fputil::multiply_add(pe, one_over_x_m, sign ? (-M_MATH_PI_2) : (M_MATH_PI_2));
  68   }
  69
  70   double pos_x = FPB(x_abs).get_val();
  71   bool one_over_x = pos_x > 1.0;
  72   if (one_over_x) {
  73     pos_x = 1.0 / pos_x;
  74   }
  75
  76   double near_x = fputil::nearest_integer(pos_x * ATAN_T_SIZE);
  77   int val = static_cast<int>(near_x);
  78   near_x *= 1.0 / ATAN_T_SIZE;
  79
  80   double v = (pos_x - near_x) / fputil::multiply_add(near_x, pos_x, 1.0);
  81   double v2 = v * v;
  82   double pe = __llvm_libc::fputil::polyeval(v2, ATAN_K[0], ATAN_K[1], ATAN_K[2],
  83                                             ATAN_K[3], ATAN_K[4]);
  84   double result;
  85   if (one_over_x)
  86     result = M_MATH_PI_2 - fputil::multiply_add(pe, v, ATAN_T[val - 1]);
  87   else
  88     result = fputil::multiply_add(pe, v, ATAN_T[val - 1]);
  89   return sign ? -result : result;
  90 }
  91
  92 // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
  93 //                 [|1, D...|], [0, 0.5]);
  94 constexpr double ASIN_COEFFS[10] = {0x1.5555555540fa1p-3, 0x1.333333512edc2p-4,
  95                                     0x1.6db6cc1541b31p-5, 0x1.f1caff324770ep-6,
  96                                     0x1.6e43899f5f4f4p-6, 0x1.1f847cf652577p-6,
  97                                     0x1.9b60f47f87146p-7, 0x1.259e2634c494fp-6,
  98                                     -0x1.df946fa875ddp-8, 0x1.02311ecf99c28p-5};
  99
 100 // Evaluate P(x^2) - 1, where P(x^2) ~ asin(x)/x
 101 LIBC_INLINE double asin_eval(double xsq) {
 102   double x4 = xsq * xsq;
 103   double r1 = fputil::polyeval(x4, ASIN_COEFFS[0], ASIN_COEFFS[2],
 104                                ASIN_COEFFS[4], ASIN_COEFFS[6], ASIN_COEFFS[8]);
 105   double r2 = fputil::polyeval(x4, ASIN_COEFFS[1], ASIN_COEFFS[3],
 106                                ASIN_COEFFS[5], ASIN_COEFFS[7], ASIN_COEFFS[9]);
 107   return fputil::multiply_add(xsq, r2, r1);
 108 }
 109
 110 } // namespace __llvm_libc
 111
 112 #endif // LLVM_LIBC_SRC_MATH_GENERIC_INV_TRIGF_UTILS_H