libc/src/math/generic/acosf.cpp

   1 //===-- Single-precision acos 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/acosf.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/sqrt.h"
  16 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
  17
  18 #include <errno.h>
  19
  20 #include "inv_trigf_utils.h"
  21
  22 namespace __llvm_libc {
  23
  24 static constexpr size_t N_EXCEPTS = 4;
  25
  26 // Exceptional values when |x| <= 0.5
  27 static constexpr fputil::ExceptValues<float, N_EXCEPTS> ACOSF_EXCEPTS = {{
  28     // (inputs, RZ output, RU offset, RD offset, RN offset)
  29     // x = 0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)
  30     {0x328885a3, 0x3fc90fda, 1, 0, 1},
  31     // x = -0x1.110b46p-26, acosf(x) = 0x1.921fb4p0 (RZ)
  32     {0xb28885a3, 0x3fc90fda, 1, 0, 1},
  33     // x = 0x1.04c444p-12, acosf(x) = 0x1.920f68p0 (RZ)
  34     {0x39826222, 0x3fc907b4, 1, 0, 1},
  35     // x = -0x1.04c444p-12, acosf(x) = 0x1.923p0 (RZ)
  36     {0xb9826222, 0x3fc91800, 1, 0, 1},
  37 }};
  38
  39 LLVM_LIBC_FUNCTION(float, acosf, (float x)) {
  40   using FPBits = typename fputil::FPBits<float>;
  41   FPBits xbits(x);
  42   uint32_t x_uint = xbits.uintval();
  43   uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
  44   uint32_t x_sign = x_uint >> 31;
  45
  46   // |x| <= 0.5
  47   if (LIBC_UNLIKELY(x_abs <= 0x3f00'0000U)) {
  48     // |x| < 0x1p-10
  49     if (LIBC_UNLIKELY(x_abs < 0x3a80'0000U)) {
  50       // When |x| < 2^-10, we use the following approximation:
  51       //   acos(x) = pi/2 - asin(x)
  52       //           ~ pi/2 - x - x^3 / 6
  53
  54       // Check for exceptional values
  55       if (auto r = ACOSF_EXCEPTS.lookup(x_uint); LIBC_UNLIKELY(r.has_value()))
  56         return r.value();
  57
  58       double xd = static_cast<double>(x);
  59       return static_cast<float>(fputil::multiply_add(
  60           -0x1.5555555555555p-3 * xd, xd * xd, M_MATH_PI_2 - xd));
  61     }
  62
  63     // For |x| <= 0.5, we approximate acosf(x) by:
  64     //   acos(x) = pi/2 - asin(x) = pi/2 - x * P(x^2)
  65     // Where P(X^2) = Q(X) is a degree-20 minimax even polynomial approximating
  66     // asin(x)/x on [0, 0.5] generated by Sollya with:
  67     // > Q = fpminimax(asin(x)/x, [|0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20|],
  68     //                 [|1, D...|], [0, 0.5]);
  69     double xd = static_cast<double>(x);
  70     double xsq = xd * xd;
  71     double x3 = xd * xsq;
  72     double r = asin_eval(xsq);
  73     return static_cast<float>(fputil::multiply_add(-x3, r, M_MATH_PI_2 - xd));
  74   }
  75
  76   // |x| > 1, return NaNs.
  77   if (LIBC_UNLIKELY(x_abs > 0x3f80'0000U)) {
  78     if (x_abs <= 0x7f80'0000U) {
  79       fputil::set_errno_if_required(EDOM);
  80       fputil::raise_except_if_required(FE_INVALID);
  81     }
  82     return x + FPBits::build_quiet_nan(0);
  83   }
  84
  85   // When 0.5 < |x| <= 1, we perform range reduction as follow:
  86   //
  87   // Assume further that 0.5 < x <= 1, and let:
  88   //   y = acos(x)
  89   // We use the double angle formula:
  90   //   x = cos(y) = 1 - 2 sin^2(y/2)
  91   // So:
  92   //   sin(y/2) = sqrt( (1 - x)/2 )
  93   // And hence:
  94   //   y = 2 * asin( sqrt( (1 - x)/2 ) )
  95   // Let u = (1 - x)/2, then
  96   //   acos(x) = 2 * asin( sqrt(u) )
  97   // Moreover, since 0.5 < x <= 1,
  98   //   0 <= u < 1/4, and 0 <= sqrt(u) < 0.5,
  99   // And hence we can reuse the same polynomial approximation of asin(x) when
 100   // |x| <= 0.5:
 101   //   acos(x) ~ 2 * sqrt(u) * P(u).
 102   //
 103   // When -1 <= x <= -0.5, we use the identity:
 104   //   acos(x) = pi - acos(-x)
 105   // which is reduced to the postive case.
 106
 107   xbits.set_sign(false);
 108   double xd = static_cast<double>(xbits.get_val());
 109   double u = fputil::multiply_add(-0.5, xd, 0.5);
 110   double cv = 2 * fputil::sqrt(u);
 111
 112   double r3 = asin_eval(u);
 113   double r = fputil::multiply_add(cv * u, r3, cv);
 114   return static_cast<float>(x_sign ? M_MATH_PI - r : r);
 115 }
 116
 117 } // namespace __llvm_libc