libc/src/math/generic/sincosf_utils.h

   1 //===-- Collection of utils for sinf/cosf/sincosf ---------------*- C++ -*-===//
   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_SINCOSF_UTILS_H
  10 #define LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF_UTILS_H
  11
  12 #include "src/__support/FPUtil/FPBits.h"
  13 #include "src/__support/FPUtil/PolyEval.h"
  14 #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
  15
  16 #if defined(LIBC_TARGET_CPU_HAS_FMA)
  17 #include "range_reduction_fma.h"
  18 // using namespace __llvm_libc::fma;
  19 using __llvm_libc::fma::FAST_PASS_BOUND;
  20 using __llvm_libc::fma::large_range_reduction;
  21 using __llvm_libc::fma::small_range_reduction;
  22 #else
  23 #include "range_reduction.h"
  24 // using namespace __llvm_libc::generic;
  25 using __llvm_libc::generic::FAST_PASS_BOUND;
  26 using __llvm_libc::generic::large_range_reduction;
  27 using __llvm_libc::generic::small_range_reduction;
  28 #endif // LIBC_TARGET_CPU_HAS_FMA
  29
  30 namespace __llvm_libc {
  31
  32 // Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
  33 // Table is generated with Sollya as follow:
  34 // > display = hexadecimal;
  35 // > for k from 0 to 63 do { D(sin(k * pi/32)); };
  36 const double SIN_K_PI_OVER_32[64] = {
  37     0x0.0000000000000p+0,  0x1.917a6bc29b42cp-4,  0x1.8f8b83c69a60bp-3,
  38     0x1.294062ed59f06p-2,  0x1.87de2a6aea963p-2,  0x1.e2b5d3806f63bp-2,
  39     0x1.1c73b39ae68c8p-1,  0x1.44cf325091dd6p-1,  0x1.6a09e667f3bcdp-1,
  40     0x1.8bc806b151741p-1,  0x1.a9b66290ea1a3p-1,  0x1.c38b2f180bdb1p-1,
  41     0x1.d906bcf328d46p-1,  0x1.e9f4156c62ddap-1,  0x1.f6297cff75cbp-1,
  42     0x1.fd88da3d12526p-1,  0x1.0000000000000p+0,  0x1.fd88da3d12526p-1,
  43     0x1.f6297cff75cbp-1,   0x1.e9f4156c62ddap-1,  0x1.d906bcf328d46p-1,
  44     0x1.c38b2f180bdb1p-1,  0x1.a9b66290ea1a3p-1,  0x1.8bc806b151741p-1,
  45     0x1.6a09e667f3bcdp-1,  0x1.44cf325091dd6p-1,  0x1.1c73b39ae68c8p-1,
  46     0x1.e2b5d3806f63bp-2,  0x1.87de2a6aea963p-2,  0x1.294062ed59f06p-2,
  47     0x1.8f8b83c69a60bp-3,  0x1.917a6bc29b42cp-4,  0x0.0000000000000p+0,
  48     -0x1.917a6bc29b42cp-4, -0x1.8f8b83c69a60bp-3, -0x1.294062ed59f06p-2,
  49     -0x1.87de2a6aea963p-2, -0x1.e2b5d3806f63bp-2, -0x1.1c73b39ae68c8p-1,
  50     -0x1.44cf325091dd6p-1, -0x1.6a09e667f3bcdp-1, -0x1.8bc806b151741p-1,
  51     -0x1.a9b66290ea1a3p-1, -0x1.c38b2f180bdb1p-1, -0x1.d906bcf328d46p-1,
  52     -0x1.e9f4156c62ddap-1, -0x1.f6297cff75cbp-1,  -0x1.fd88da3d12526p-1,
  53     -0x1.0000000000000p+0, -0x1.fd88da3d12526p-1, -0x1.f6297cff75cbp-1,
  54     -0x1.e9f4156c62ddap-1, -0x1.d906bcf328d46p-1, -0x1.c38b2f180bdb1p-1,
  55     -0x1.a9b66290ea1a3p-1, -0x1.8bc806b151741p-1, -0x1.6a09e667f3bcdp-1,
  56     -0x1.44cf325091dd6p-1, -0x1.1c73b39ae68c8p-1, -0x1.e2b5d3806f63bp-2,
  57     -0x1.87de2a6aea963p-2, -0x1.294062ed59f06p-2, -0x1.8f8b83c69a60bp-3,
  58     -0x1.917a6bc29b42cp-4,
  59 };
  60
  61 LIBC_INLINE void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,
  62                               double &cos_k, double &sin_y, double &cosm1_y) {
  63   int64_t k;
  64   double y;
  65
  66   if (LIBC_LIKELY(x_abs < FAST_PASS_BOUND)) {
  67     k = small_range_reduction(xd, y);
  68   } else {
  69     fputil::FPBits<float> x_bits(x_abs);
  70     k = large_range_reduction(xd, x_bits.get_exponent(), y);
  71   }
  72
  73   // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
  74   // So k is an integer and -0.5 <= y <= 0.5.
  75   // Then sin(x) = sin((k + y)*pi/32)
  76   //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
  77
  78   sin_k = SIN_K_PI_OVER_32[k & 63];
  79   // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
  80   // cos_k = cos(k * pi/32)
  81   cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
  82
  83   double ysq = y * y;
  84
  85   // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
  86   // with:
  87   // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
  88   sin_y =
  89       y * fputil::polyeval(ysq, 0x1.921fb54442d18p-4, -0x1.4abbce625abb1p-13,
  90                            0x1.466bc624f2776p-24, -0x1.32c3a619d4a7ep-36);
  91   // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
  92   // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
  93   // Note that cosm1_y = cos(y*pi/32) - 1.
  94   cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3cc9be430bp-8,
  95                                    0x1.03c1f070c2e27p-18, -0x1.55cc84bd942p-30);
  96 }
  97
  98 } // namespace __llvm_libc
  99
 100 #endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF_UTILS_H