[libc] Change sinf/cosf range reduction to mod pi/32 to be shared with tanf.

[lntue]

Summary

Change sinf/cosf range reduction to mod pi/32 to be shared with tanf,
since polynomial approximations for tanf on subintervals of length pi/16 do not
provide enough accuracy.

[lntue]

Small fixes for comments.

libc/docs/math.rst

* We also use the *perf* tool from the `CORE-MATH <https://core-math.gitlabpages.inria.fr/>`_
   table below, using the system library as reference (such as the `GNU C library <https://www.gnu.org/software/libc/>`_
   on Linux). Fmod performance results obtained with "differential_testing".
 
 +--------------+-------------------------------+-------------------------------+-------------------------------------+---------------------------------------------------------------------+
 | <Func>       | Reciprocal throughput (ns)    | Latency (ns)                  | Testing ranges                      | Testing configuration                                               |
 |              +-----------+-------------------+-----------+-------------------+                                     +------------+-------------------------+--------------+---------------+
 |              | LLVM libc | Reference (glibc) | LLVM libc | Reference (glibc) |                                     | CPU        | OS                      | Compiler     | Special flags |
 +==============+===========+===================+===========+===================+=====================================+============+=========================+==============+===============+
-| cosf         |        14 |                32 |        56 |                59 | :math:`[0, 2\pi]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
+| cosf         |        13 |                32 |        53 |                59 | :math:`[0, 2\pi]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | coshf        |        23 |                20 |        73 |                49 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | expf         |         9 |                 7 |        44 |                38 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | exp2f        |         8 |                 6 |        35 |                23 | :math:`[-10, 10]`                   | i5-1135G7  | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +              +-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 |              |        11 |                 6 |        49 |                31 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |

 | logf         |        12 |                10 |        56 |                46 | :math:`[e^{-1}, e]`                 | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | log10f       |        13 |                25 |        57 |                72 | :math:`[e^{-1}, e]`                 | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | log1pf       |        16 |                33 |        61 |                97 | :math:`[e^{-0.5} - 1, e^{0.5} - 1]` | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | log2f        |        13 |                10 |        57 |                46 | :math:`[e^{-1}, e]`                 | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
-| sinf         |        13 |                25 |        54 |                57 | :math:`[-\pi, \pi]`                 | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
+| sinf         |        12 |                25 |        51 |                57 | :math:`[-\pi, \pi]`                 | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
-| sincosf      |        20 |                30 |        62 |                68 | :math:`[-\pi, \pi]`                 | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
+| sincosf      |        19 |                30 |        57 |                68 | :math:`[-\pi, \pi]`                 | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | sinhf        |        23 |                64 |        73 |               141 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | tanhf        |        25 |                59 |        95 |               125 | :math:`[-10, 10]`                   | Ryzen 1700 | Ubuntu 20.04 LTS x86_64 | Clang 12.0.0 | FMA           |
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 
 References
 ==========
 
 * `CRLIBM <https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804/file/crlibm.pdf>`_.
 * `RLIBM <https://people.cs.rutgers.edu/~sn349/rlibm/>`_.
 * `Sollya <https://www.sollya.org/>`_.
 * `The CORE-MATH Project <https://core-math.gitlabpages.inria.fr/>`_.
 * `The GNU C Library (glibc) <https://www.gnu.org/software/libc/>`_.
 * `The GNU MPFR Library <https://www.mpfr.org/>`_.

libc/src/math/generic/cosf.cpp

LLVM_LIBC_FUNCTION(float, cosf, (float x)) {
   xbits.set_sign(false);
 
   uint32_t x_abs = xbits.uintval();
   double xd = static_cast<double>(xbits.get_val());
 
   // Range reduction:
   // For |x| > pi/16, we perform range reduction as follows:
   // Find k and y such that:
-  //   x = (k + y) * pi/16
+  //   x = (k + y) * pi/32
   //   k is an integer
   //   |y| < 0.5
-  // For small range (|x| < 2^46 when FMA instructions are available, 2^22
+  // For small range (|x| < 2^45 when FMA instructions are available, 2^22
   // otherwise), this is done by performing:
-  //   k = round(x * 16/pi)
-  //   y = x * 16/pi - k
+  //   k = round(x * 32/pi)
+  //   y = x * 32/pi - k
   // For large range, we will omit all the higher parts of 16/pi such that the
-  // least significant bits of their full products with x are larger than 31,
-  // since cos((k + y + 32*i) * pi/16) = cos(x + i * 2pi) = cos(x).
+  // least significant bits of their full products with x are larger than 63,
+  // since cos((k + y + 64*i) * pi/32) = cos(x + i * 2pi) = cos(x).
   //
-  // When FMA instructions are not available, we store the digits of 16/pi in
+  // When FMA instructions are not available, we store the digits of 32/pi in
   // chunks of 28-bit precision.  This will make sure that the products:
-  //   x * SIXTEEN_OVER_PI_28[i] are all exact.
-  // When FMA instructions are available, we simply store the digits of 16/pi in
+  //   x * THIRTYTWO_OVER_PI_28[i] are all exact.
+  // When FMA instructions are available, we simply store the digits of 32/pi in
   // chunks of doubles (53-bit of precision).
   // So when multiplying by the largest values of single precision, the
   // resulting output should be correct up to 2^(-208 + 128) ~ 2^-80.  By the
   // worst-case analysis of range reduction, |y| >= 2^-38, so this should give
   // us more than 40 bits of accuracy. For the worst-case estimation of range
   // reduction, see for instances:
   //   Elementary Functions by J-M. Muller, Chapter 11,
   //   Handbook of Floating-Point Arithmetic by J-M. Muller et. al.,
   //   Chapter 10.2.
   //
   // Once k and y are computed, we then deduce the answer by the cosine of sum
   // formula:
-  //   cos(x) = cos((k + y)*pi/16)
-  //          = cos(y*pi/16) * cos(k*pi/16) - sin(y*pi/16) * sin(k*pi/16)
-  // The values of sin(k*pi/16) and cos(k*pi/16) for k = 0..31 are precomputed
-  // and stored using a vector of 32 doubles. Sin(y*pi/16) and cos(y*pi/16) are
-  // computed using degree-7 and degree-8 minimax polynomials generated by
+  //   cos(x) = cos((k + y)*pi/32)
+  //          = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
+  // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..63 are precomputed
+  // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
+  // computed using degree-7 and degree-6 minimax polynomials generated by
   // Sollya respectively.
 
   // |x| < 0x1.0p-12f
   if (unlikely(x_abs < 0x3980'0000U)) {
     // When |x| < 2^-12, the relative error of the approximation cos(x) ~ 1
     // is:
     //   |cos(x) - 1| < |x^2 / 2| = 2^-25 < epsilon(1)/2.
     // So the correctly rounded values of cos(x) are:
if (x_abs == 0x7f80'0000U) {
       errno = EDOM;
       fputil::set_except(FE_INVALID);
     }
     return x +
            FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
   }
 
   // Combine the results with the sine of sum formula:
-  //   cos(x) = cos((k + y)*pi/16)
-  //          = cos(y*pi/16) * cos(k*pi/16) - sin(y*pi/16) * sin(k*pi/16)
+  //   cos(x) = cos((k + y)*pi/32)
+  //          = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
   //          = cosm1_y * cos_k + sin_y * sin_k
   //          = (cosm1_y * cos_k + cos_k) + sin_y * sin_k
   double sin_k, cos_k, sin_y, cosm1_y;
 
   sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
 
   return fputil::multiply_add(sin_y, -sin_k,
                               fputil::multiply_add(cosm1_y, cos_k, cos_k));
 }
 
 } // namespace __llvm_libc

libc/src/math/generic/range_reduction.h

 namespace __llvm_libc {
 
 namespace generic {
 
 static constexpr uint32_t FAST_PASS_BOUND = 0x4a80'0000U; // 2^22
 
 static constexpr int N_ENTRIES = 8;
 
-// We choose to split bits of 16/pi into 28-bit precision pieces, so that the
-// product of x * SIXTEEN_OVER_PI_28[i] is exact.
+// We choose to split bits of 32/pi into 28-bit precision pieces, so that the
+// product of x * THIRTYTWO_OVER_PI_28[i] is exact.
 // These are generated by Sollya with:
-// > a1 = D(round(16/pi, 28, RN)); a1;
-// > a2 = D(round(16/pi - a1, 28, RN)); a2;
-// > a3 = D(round(16/pi - a1 - a2, 28, RN)); a3;
-// > a4 = D(round(16/pi - a1 - a2 - a3, 28, RN)); a4;
+// > a1 = D(round(32/pi, 28, RN)); a1;
+// > a2 = D(round(32/pi - a1, 28, RN)); a2;
+// > a3 = D(round(32/pi - a1 - a2, 28, RN)); a3;
+// > a4 = D(round(32/pi - a1 - a2 - a3, 28, RN)); a4;
 // ...
-static constexpr double SIXTEEN_OVER_PI_28[N_ENTRIES] = {
-    0x1.45f306ep+2,   -0x1.b1bbeaep-29,  0x1.3f84ebp-58,    -0x1.7056592p-88,
-    0x1.c0db62ap-117, -0x1.4cd8778p-146, -0x1.bef806cp-175, 0x1.63abdecp-205};
+static constexpr double THIRTYTWO_OVER_PI_28[N_ENTRIES] = {
+    0x1.45f306ep+3,   -0x1.b1bbeaep-28,  0x1.3f84ebp-57,    -0x1.7056592p-87,
+    0x1.c0db62ap-116, -0x1.4cd8778p-145, -0x1.bef806cp-174, 0x1.63abdecp-204};
 
 // Exponents of the least significant bits of the corresponding entries in
-// SIXTEEN_OVER_PI_28.
-static constexpr int SIXTEEN_OVER_PI_28_LSB_EXP[N_ENTRIES] = {
-    -25, -56, -82, -115, -144, -171, -201, -231};
+// THIRTYTWO_OVER_PI_28.
+static constexpr int THIRTYTWO_OVER_PI_28_LSB_EXP[N_ENTRIES] = {
+    -24, -55, -81, -114, -143, -170, -200, -230};
 
 // Return k and y, where
 //   k = round(x * 16 / pi) and y = (x * 16 / pi) - k.
 static inline int64_t small_range_reduction(double x, double &y) {
-  double prod = x * SIXTEEN_OVER_PI_28[0];
+  double prod = x * THIRTYTWO_OVER_PI_28[0];
   double kd = fputil::nearest_integer(prod);
   y = prod - kd;
-  y = fputil::multiply_add(x, SIXTEEN_OVER_PI_28[1], y);
-  y = fputil::multiply_add(x, SIXTEEN_OVER_PI_28[2], y);
+  y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[1], y);
+  y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[2], y);
   return static_cast<int64_t>(kd);
 }
 
 // Return k and y, where
-//   k = round(x * 16 / pi) and y = (x * 16 / pi) - k.
-// For large range, there are at most 2 parts of SIXTEEN_OVER_PI_28 contributing
-// to the lowest 5 binary digits (k & 31).  If the least significant bit of
-// x * the least significant bit of SIXTEEN_OVER_PI_28[i] >= 32, we can
-// completely ignore SIXTEEN_OVER_PI_28[i].
+//   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
+// For large range, there are at most 2 parts of THIRTYTWO_OVER_PI_28
+// contributing to the lowest 6 binary digits (k & 63).  If the least
+// significant bit of x * the least significant bit of THIRTYTWO_OVER_PI_28[i]
+// >= 64, we can completely ignore THIRTYTWO_OVER_PI_28[i].
 static inline int64_t large_range_reduction(double x, int x_exp, double &y) {
   int idx = 0;
   y = 0;
   int x_lsb_exp_m4 = x_exp - fputil::FloatProperties<float>::MANTISSA_WIDTH;
 
-  // Skipping the first parts of 16/pi such that:
-  //   LSB of x * LSB of SIXTEEN_OVER_PI_28[i] >= 32.
-  while (x_lsb_exp_m4 + SIXTEEN_OVER_PI_28_LSB_EXP[idx] > 4)
+  // Skipping the first parts of 32/pi such that:
+  //   LSB of x * LSB of THIRTYTWO_OVER_PI_28[i] >= 32.
+  while (x_lsb_exp_m4 + THIRTYTWO_OVER_PI_28_LSB_EXP[idx] > 5)
     ++idx;
 
-  double prod_hi = x * SIXTEEN_OVER_PI_28[idx];
-  // Get the integral part of x * SIXTEEN_OVER_PI_28[idx]
+  double prod_hi = x * THIRTYTWO_OVER_PI_28[idx];
+  // Get the integral part of x * THIRTYTWO_OVER_PI_28[idx]
   double k_hi = fputil::nearest_integer(prod_hi);
-  // Get the fractional part of x * SIXTEEN_OVER_PI_28[idx]
+  // Get the fractional part of x * THIRTYTWO_OVER_PI_28[idx]
   double frac = prod_hi - k_hi;
-  double prod_lo = fputil::multiply_add(x, SIXTEEN_OVER_PI_28[idx + 1], frac);
+  double prod_lo = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 1], frac);
   double k_lo = fputil::nearest_integer(prod_lo);
 
   // Now y is the fractional parts.
   y = prod_lo - k_lo;
-  y = fputil::multiply_add(x, SIXTEEN_OVER_PI_28[idx + 2], y);
-  y = fputil::multiply_add(x, SIXTEEN_OVER_PI_28[idx + 3], y);
+  y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 2], y);
+  y = fputil::multiply_add(x, THIRTYTWO_OVER_PI_28[idx + 3], y);
 
   return static_cast<int64_t>(k_hi) + static_cast<int64_t>(k_lo);
 }
 
-// Exceptional cases.
-static constexpr int N_EXCEPTS = 3;
-
-static constexpr fputil::ExceptionalValues<float, N_EXCEPTS> SinfExcepts{
-    /* inputs */ {
-        0x3fa7832a, // x = 0x1.4f0654p0
-        0x46199998, // x = 0x1.33333p13
-        0x55cafb2a, // x = 0x1.95f654p44
-    },
-    /* outputs (RZ, RU offset, RD offset, RN offset) */
-    {
-        {0x3f7741b5, 1, 0, 1}, // x = 0x1.4f0654p0, sin(x) = 0x1.ee836ap-1 (RZ)
-        {0xbeb1fa5d, 0, 1, 0}, // x = 0x1.33333p13, sin(x) = -0x1.63f4bap-2 (RZ)
-        {0xbf7e7a16, 0, 1,
-         1}, // x = 0x1.95f654p44, sin(x) = -0x1.fcf42cp-1 (RZ)
-    }};
-
 } // namespace generic
 
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_H

libc/src/math/generic/range_reduction_fma.h

 //===-- Utilities for trigonometric functions with FMA ----------*- C++ -*-===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H
 #define LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H
 
 #include "src/__support/FPUtil/FMA.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/except_value_utils.h"
 #include "src/__support/FPUtil/nearest_integer.h"
 
+// #include <iostream>
+
 namespace __llvm_libc {
 
 namespace fma {
 
-static constexpr uint32_t FAST_PASS_BOUND = 0x5680'0000U; // 2^46
+static constexpr uint32_t FAST_PASS_BOUND = 0x5600'0000U; // 2^45
 
-// Digits of 1/pi, generated by Sollya with:
-// > a0 = D(16/pi);
-// > a1 = D(16/pi - a0);
-// > a2 = D(16/pi - a0 - a1);
-// > a3 = D(16/pi - a0 - a1 - a2);
-static constexpr double SIXTEEN_OVER_PI[5] = {
-    0x1.45f306dc9c883p+2, -0x1.6b01ec5417056p-52, -0x1.6447e493ad4cep-106,
-    0x1.e21c820ff28b2p-160, -0x1.508510ea79237p-215};
+// Digits of 32/pi, generated by Sollya with:
+// > a0 = D(32/pi);
+// > a1 = D(32/pi - a0);
+// > a2 = D(32/pi - a0 - a1);
+// > a3 = D(32/pi - a0 - a1 - a2);
+static constexpr double THIRTYTWO_OVER_PI[5] = {
+    0x1.45f306dc9c883p+3, -0x1.6b01ec5417056p-51, -0x1.6447e493ad4cep-105,
+    0x1.e21c820ff28b2p-159, -0x1.508510ea79237p-214};
 
 // Return k and y, where
-//   k = round(x * 16 / pi) and y = (x * 16 / pi) - k.
-// Assume x is non-negative.
+//   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
 static inline int64_t small_range_reduction(double x, double &y) {
-  double kd = fputil::nearest_integer(x * SIXTEEN_OVER_PI[0]);
-  y = fputil::fma(x, SIXTEEN_OVER_PI[0], -kd);
-  y = fputil::fma(x, SIXTEEN_OVER_PI[1], y);
+  // std::cout << "  We should NOT be here!!!" << std::endl;
+  double kd = fputil::nearest_integer(x * THIRTYTWO_OVER_PI[0]);
+  y = fputil::fma(x, THIRTYTWO_OVER_PI[0], -kd);
+  y = fputil::fma(x, THIRTYTWO_OVER_PI[1], y);
   return static_cast<int64_t>(kd);
 }
 
 // Return k and y, where
-//   k = round(x * 16 / pi) and y = (x * 16 / pi) - k.
+//   k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
+// This is used for sinf, cosf, sincosf.
 static inline int64_t large_range_reduction(double x, int x_exp, double &y) {
-  // 2^46 <= |x| < 2^99
-  if (x_exp < 99) {
-    // - When x < 2^99, the full exact product of x * SIXTEEN_OVER_PI[0]
-    // contains at least one integral bit <= 2^4.
-    // - When 2^46 <= |x| < 2^56, the lowest 5 unit bits are contained
-    // in the last 10 bits of double(x * SIXTEEN_OVER_PI[0]).
-    // - When |x| >= 2^56, the LSB of double(x * SIXTEEN_OVER_PI[0]) is at least
-    // 32.
-    fputil::FPBits<double> prod_hi(x * SIXTEEN_OVER_PI[0]);
-    prod_hi.bits &= (x_exp < 56) ? (~0xfffULL) : (~0ULL); // |x| < 2^56
+  // 2^45 <= |x| < 2^100
+  if (x_exp < 100) {
+    // - When x < 2^100, the full exact product of x * THIRTYTWO_OVER_PI[0]
+    // contains at least one integral bit <= 2^5.
+    // - When 2^45 <= |x| < 2^55, the lowest 6 unit bits are contained
+    // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[0]).
+    // - When |x| >= 2^55, the LSB of double(x * THIRTYTWO_OVER_PI[0]) is at
+    // least 2^6.
+    fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[0]);
+    prod_hi.bits &= (x_exp < 55) ? (~0xfffULL) : (~0ULL); // |x| < 2^55
     double k_hi = fputil::nearest_integer(static_cast<double>(prod_hi));
-    double truncated_prod = fputil::fma(x, SIXTEEN_OVER_PI[0], -k_hi);
-    double prod_lo = fputil::fma(x, SIXTEEN_OVER_PI[1], truncated_prod);
+    double truncated_prod = fputil::fma(x, THIRTYTWO_OVER_PI[0], -k_hi);
+    double prod_lo = fputil::fma(x, THIRTYTWO_OVER_PI[1], truncated_prod);
     double k_lo = fputil::nearest_integer(prod_lo);
-    y = fputil::fma(x, SIXTEEN_OVER_PI[1], truncated_prod - k_lo);
-    y = fputil::fma(x, SIXTEEN_OVER_PI[2], y);
-    y = fputil::fma(x, SIXTEEN_OVER_PI[3], y);
+    y = fputil::fma(x, THIRTYTWO_OVER_PI[1], truncated_prod - k_lo);
+    y = fputil::fma(x, THIRTYTWO_OVER_PI[2], y);
+    y = fputil::fma(x, THIRTYTWO_OVER_PI[3], y);
 
     return static_cast<int64_t>(k_lo);
   }
 
-  // - When x >= 2^110, the full exact product of x * SIXTEEN_OVER_PI[0] does
-  // not contain any of the lowest 5 unit bits, so we can ignore it completely.
-  // - When 2^99 <= |x| < 2^110, the lowest 5 unit bits are contained
-  // in the last 12 bits of double(x * SIXTEEN_OVER_PI[1]).
-  // - When |x| >= 2^110, the LSB of double(x * SIXTEEN_OVER_PI[1]) is at
-  // least 32.
-  fputil::FPBits<double> prod_hi(x * SIXTEEN_OVER_PI[1]);
+  // - When x >= 2^110, the full exact product of x * THIRTYTWO_OVER_PI[0] does
+  // not contain any of the lowest 6 unit bits, so we can ignore it completely.
+  // - When 2^100 <= |x| < 2^110, the lowest 6 unit bits are contained
+  // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[1]).
+  // - When |x| >= 2^110, the LSB of double(x * THIRTYTWO_OVER_PI[1]) is at
+  // least 64.
+  fputil::FPBits<double> prod_hi(x * THIRTYTWO_OVER_PI[1]);
   prod_hi.bits &= (x_exp < 110) ? (~0xfffULL) : (~0ULL); // |x| < 2^110
   double k_hi = fputil::nearest_integer(static_cast<double>(prod_hi));
-  double truncated_prod = fputil::fma(x, SIXTEEN_OVER_PI[1], -k_hi);
-  double prod_lo = fputil::fma(x, SIXTEEN_OVER_PI[2], truncated_prod);
+  double truncated_prod = fputil::fma(x, THIRTYTWO_OVER_PI[1], -k_hi);
+  double prod_lo = fputil::fma(x, THIRTYTWO_OVER_PI[2], truncated_prod);
   double k_lo = fputil::nearest_integer(prod_lo);
-  y = fputil::fma(x, SIXTEEN_OVER_PI[2], truncated_prod - k_lo);
-  y = fputil::fma(x, SIXTEEN_OVER_PI[3], y);
-  y = fputil::fma(x, SIXTEEN_OVER_PI[4], y);
+  y = fputil::fma(x, THIRTYTWO_OVER_PI[2], truncated_prod - k_lo);
+  y = fputil::fma(x, THIRTYTWO_OVER_PI[3], y);
+  y = fputil::fma(x, THIRTYTWO_OVER_PI[4], y);
 
   return static_cast<int64_t>(k_lo);
 }
 
-// Exceptional cases.
-static constexpr int N_EXCEPTS = 2;
-
-static constexpr fputil::ExceptionalValues<float, N_EXCEPTS> SinfExcepts{
-    /* inputs */ {
-        0x3fa7832a, // x = 0x1.4f0654p0
-        0x55cafb2a, // x = 0x1.95f654p44
-    },
-    /* outputs (RZ, RU offset, RD offset, RN offset) */
-    {
-        {0x3f7741b5, 1, 0, 1}, // x = 0x1.4f0654p0, sin(x) = 0x1.ee836ap-1 (RZ)
-        {0xbf7e7a16, 0, 1,
-         1}, // x = 0x1.95f654p44, sin(x) = -0x1.fcf42cp-1 (RZ)
-    }};
-
 } // namespace fma
 
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_RANGE_REDUCTION_FMA_H

libc/src/math/generic/sincosf.cpp

 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/common.h"
 
 #include <errno.h>
 
 namespace __llvm_libc {
 
 // Exceptional values
-static constexpr int N_EXCEPTS = 10;
+static constexpr int N_EXCEPTS = 6;
 
 static constexpr uint32_t EXCEPT_INPUTS[N_EXCEPTS] = {
-    0x3b5637f5, // x = 0x1.ac6feap-9
-    0x3fa7832a, // x = 0x1.4f0654p0
-    0x46199998, // x = 0x1.33333p13
-    0x55325019, // x = 0x1.64a032p43
-    0x55cafb2a, // x = 0x1.95f654p44
-    0x5922aa80, // x = 0x1.4555p51
-    0x5aa4542c, // x = 0x1.48a858p54
-    0x5f18b878, // x = 0x1.3170fp63
-    0x6115cb11, // x = 0x1.2b9622p67
-    0x7beef5ef, // x = 0x1.ddebdep120
+    0x46199998, // x = 0x1.33333p13   x
+    0x55325019, // x = 0x1.64a032p43  x
+    0x5922aa80, // x = 0x1.4555p51    x
+    0x5f18b878, // x = 0x1.3170fp63   x
+    0x6115cb11, // x = 0x1.2b9622p67  x
+    0x7beef5ef, // x = 0x1.ddebdep120 x
 };
 
 static constexpr uint32_t EXCEPT_OUTPUTS_SIN[N_EXCEPTS][4] = {
-    {0x3b5637dc, 1, 0, 0}, // x = 0x1.ac6feap-9, sin(x) = 0x1.ac6fb8p-9 (RZ)
-    {0x3f7741b5, 1, 0, 1}, // x = 0x1.4f0654p0, sin(x) = 0x1.ee836ap-1 (RZ)
     {0xbeb1fa5d, 0, 1, 0}, // x = 0x1.33333p13, sin(x) = -0x1.63f4bap-2 (RZ)
     {0xbf171adf, 0, 1, 1}, // x = 0x1.64a032p43, sin(x) = -0x1.2e35bep-1 (RZ)
-    {0xbf7e7a16, 0, 1, 1}, // x = 0x1.95f654p44, sin(x) = -0x1.fcf42cp-1 (RZ)
     {0xbf587521, 0, 1, 1}, // x = 0x1.4555p51, sin(x) = -0x1.b0ea42p-1 (RZ)
-    {0x3f5f5646, 1, 0, 0}, // x = 0x1.48a858p54, sin(x) = 0x1.beac8cp-1 (RZ)
     {0x3dad60f6, 1, 0, 1}, // x = 0x1.3170fp63, sin(x) = 0x1.5ac1ecp-4 (RZ)
     {0xbe7cc1e0, 0, 1, 1}, // x = 0x1.2b9622p67, sin(x) = -0x1.f983cp-3 (RZ)
     {0xbf587d1b, 0, 1, 1}, // x = 0x1.ddebdep120, sin(x) = -0x1.b0fa36p-1 (RZ)
 };
 
 static constexpr uint32_t EXCEPT_OUTPUTS_COS[N_EXCEPTS][4] = {
-    {0x3f7fffa6, 1, 0, 0}, // x = 0x1.ac6feap-9, cos(x) = 0x1.ffff4cp-1 (RZ)
-    {0x3e84aabf, 1, 0, 1}, // x = 0x1.4f0654p0, cos(x) = 0x1.09557ep-2 (RZ)
     {0xbf70090b, 0, 1, 0}, // x = 0x1.33333p13, cos(x) = -0x1.e01216p-1 (RZ)
     {0x3f4ea5d2, 1, 0, 0}, // x = 0x1.64a032p43, cos(x) = 0x1.9d4ba4p-1 (RZ)
-    {0x3ddf11f3, 1, 0, 1}, // x = 0x1.95f654p44, cos(x) = 0x1.be23e6p-4 (RZ)
     {0x3f08aebe, 1, 0, 1}, // x = 0x1.4555p51, cos(x) = 0x1.115d7cp-1 (RZ)
-    {0x3efa40a4, 1, 0, 0}, // x = 0x1.48a858p54, cos(x) = 0x1.f48148p-2 (RZ)
     {0x3f7f14bb, 1, 0, 0}, // x = 0x1.3170fp63, cos(x) = 0x1.fe2976p-1 (RZ)
     {0x3f78142e, 1, 0, 1}, // x = 0x1.2b9622p67, cos(x) = 0x1.f0285cp-1 (RZ)
     {0x3f08a21c, 1, 0, 0}, // x = 0x1.ddebdep120, cos(x) = 0x1.114438p-1 (RZ)
 };
 
-// Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative
-// error is 0.5303 * 2^-23. A single-step range reduction is used for
-// small values. Large inputs have their range reduced using fast integer
-// arithmetic.
 LLVM_LIBC_FUNCTION(void, sincosf, (float x, float *sinp, float *cosp)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
 
   uint32_t x_abs = xbits.uintval() & 0x7fff'ffffU;
   double xd = static_cast<double>(x);
 
   // Range reduction:
-  // For |x| > pi/16, we perform range reduction as follows:
+  // For |x| >= 2^-12, we perform range reduction as follows:
   // Find k and y such that:
-  //   x = (k + y) * pi/16
+  //   x = (k + y) * pi/32
   //   k is an integer
   //   |y| < 0.5
   // For small range (|x| < 2^46 when FMA instructions are available, 2^22
   // otherwise), this is done by performing:
-  //   k = round(x * 16/pi)
-  //   y = x * 16/pi - k
-  // For large range, we will omit all the higher parts of 16/pi such that the
-  // least significant bits of their full products with x are larger than 31,
+  //   k = round(x * 32/pi)
+  //   y = x * 32/pi - k
+  // For large range, we will omit all the higher parts of 32/pi such that the
+  // least significant bits of their full products with x are larger than 63,
   // since:
-  //     sin((k + y + 32*i) * pi/16) = sin(x + i * 2pi) = sin(x), and
-  //     cos((k + y + 32*i) * pi/16) = cos(x + i * 2pi) = cos(x).
+  //     sin((k + y + 64*i) * pi/32) = sin(x + i * 2pi) = sin(x), and
+  //     cos((k + y + 64*i) * pi/32) = cos(x + i * 2pi) = cos(x).
   //
-  // When FMA instructions are not available, we store the digits of 16/pi in
+  // When FMA instructions are not available, we store the digits of 32/pi in
   // chunks of 28-bit precision.  This will make sure that the products:
-  //   x * SIXTEEN_OVER_PI_28[i] are all exact.
-  // When FMA instructions are available, we simply store the digits of 16/pi in
+  //   x * THIRTYTWO_OVER_PI_28[i] are all exact.
+  // When FMA instructions are available, we simply store the digits of326/pi in
   // chunks of doubles (53-bit of precision).
   // So when multiplying by the largest values of single precision, the
   // resulting output should be correct up to 2^(-208 + 128) ~ 2^-80.  By the
   // worst-case analysis of range reduction, |y| >= 2^-38, so this should give
   // us more than 40 bits of accuracy. For the worst-case estimation of range
   // reduction, see for instances:
   //   Elementary Functions by J-M. Muller, Chapter 11,
   //   Handbook of Floating-Point Arithmetic by J-M. Muller et. al.,
   //   Chapter 10.2.
   //
   // Once k and y are computed, we then deduce the answer by the sine and cosine
   // of sum formulas:
-  //   sin(x) = sin((k + y)*pi/16)
-  //          = sin(y*pi/16) * cos(k*pi/16) + cos(y*pi/16) * sin(k*pi/16)
-  //   cos(x) = cos((k + y)*pi/16)
-  //          = cos(y*pi/16) * cos(k*pi/16) - sin(y*pi/16) * sin(k*pi/16)
-  // The values of sin(k*pi/16) and cos(k*pi/16) for k = 0..31 are precomputed
-  // and stored using a vector of 32 doubles. Sin(y*pi/16) and cos(y*pi/16) are
-  // computed using degree-7 and degree-8 minimax polynomials generated by
+  //   sin(x) = sin((k + y)*pi/32)
+  //          = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
+  //   cos(x) = cos((k + y)*pi/32)
+  //          = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
+  // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..63 are precomputed
+  // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
+  // computed using degree-7 and degree-6 minimax polynomials generated by
   // Sollya respectively.
 
   // |x| < 0x1.0p-12f
   if (unlikely(x_abs < 0x3980'0000U)) {
     if (unlikely(x_abs == 0U)) {
       // For signed zeros.
       *sinp = x;
       *cosp = 1.0f;
if (unlikely(x_abs == EXCEPT_INPUTS[i])) {
       *sinp = x_sign ? -FPBits(s).get_val() : FPBits(s).get_val();
       *cosp = FPBits(c).get_val();
 
       return;
     }
   }
 
   // Combine the results with the sine and cosine of sum formulas:
-  //   sin(x) = sin((k + y)*pi/16)
-  //          = sin(y*pi/16) * cos(k*pi/16) + cos(y*pi/16) * sin(k*pi/16)
+  //   sin(x) = sin((k + y)*pi/32)
+  //          = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
   //          = sin_y * cos_k + (1 + cosm1_y) * sin_k
   //          = sin_y * cos_k + (cosm1_y * sin_k + sin_k)
-  //   cos(x) = cos((k + y)*pi/16)
-  //          = cos(y*pi/16) * cos(k*pi/16) - sin(y*pi/16) * sin(k*pi/16)
+  //   cos(x) = cos((k + y)*pi/32)
+  //          = cos(y*pi/32) * cos(k*pi/32) - sin(y*pi/32) * sin(k*pi/32)
   //          = cosm1_y * cos_k + sin_y * sin_k
   //          = (cosm1_y * cos_k + cos_k) + sin_y * sin_k
   double sin_k, cos_k, sin_y, cosm1_y;
 
   sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
 
   *sinp = fputil::multiply_add(sin_y, cos_k,
                                fputil::multiply_add(cosm1_y, sin_k, sin_k));
   *cosp = fputil::multiply_add(sin_y, -sin_k,
                                fputil::multiply_add(cosm1_y, cos_k, cos_k));
 }
 
 } // namespace __llvm_libc

libc/src/math/generic/sincosf_utils.h

 // using namespace __llvm_libc::generic;
 using __llvm_libc::generic::FAST_PASS_BOUND;
 using __llvm_libc::generic::large_range_reduction;
 using __llvm_libc::generic::small_range_reduction;
 #endif // LIBC_TARGET_HAS_FMA
 
 namespace __llvm_libc {
 
-// Lookup table for sin(k * pi / 16) with k = 0, ..., 31.
+// Lookup table for sin(k * pi / 32) with k = 0, ..., 63.
 // Table is generated with Sollya as follow:
 // > display = hexadecimal;
-// > for k from 0 to 31 do { D(sin(k * pi/16)); };
-const double SIN_K_PI_OVER_16[32] = {
-    0x0.0000000000000p+0,  0x1.8f8b83c69a60bp-3,  0x1.87de2a6aea963p-2,
-    0x1.1c73b39ae68c8p-1,  0x1.6a09e667f3bcdp-1,  0x1.a9b66290ea1a3p-1,
-    0x1.d906bcf328d46p-1,  0x1.f6297cff75cb0p-1,  0x1.0000000000000p+0,
-    0x1.f6297cff75cb0p-1,  0x1.d906bcf328d46p-1,  0x1.a9b66290ea1a3p-1,
-    0x1.6a09e667f3bcdp-1,  0x1.1c73b39ae68c8p-1,  0x1.87de2a6aea963p-2,
-    0x1.8f8b83c69a60bp-3,  0x0.0000000000000p+0,  -0x1.8f8b83c69a60bp-3,
-    -0x1.87de2a6aea963p-2, -0x1.1c73b39ae68c8p-1, -0x1.6a09e667f3bcdp-1,
-    -0x1.a9b66290ea1a3p-1, -0x1.d906bcf328d46p-1, -0x1.f6297cff75cb0p-1,
-    -0x1.0000000000000p+0, -0x1.f6297cff75cb0p-1, -0x1.d906bcf328d46p-1,
-    -0x1.a9b66290ea1a3p-1, -0x1.6a09e667f3bcdp-1, -0x1.1c73b39ae68c8p-1,
-    -0x1.87de2a6aea963p-2, -0x1.8f8b83c69a60bp-3};
+// > for k from 0 to 63 do { D(sin(k * pi/32)); };
+const double SIN_K_PI_OVER_32[64] = {
+    0x0.0000000000000p+0,  0x1.917a6bc29b42cp-4,  0x1.8f8b83c69a60bp-3,
+    0x1.294062ed59f06p-2,  0x1.87de2a6aea963p-2,  0x1.e2b5d3806f63bp-2,
+    0x1.1c73b39ae68c8p-1,  0x1.44cf325091dd6p-1,  0x1.6a09e667f3bcdp-1,
+    0x1.8bc806b151741p-1,  0x1.a9b66290ea1a3p-1,  0x1.c38b2f180bdb1p-1,
+    0x1.d906bcf328d46p-1,  0x1.e9f4156c62ddap-1,  0x1.f6297cff75cbp-1,
+    0x1.fd88da3d12526p-1,  0x1.0000000000000p+0,  0x1.fd88da3d12526p-1,
+    0x1.f6297cff75cbp-1,   0x1.e9f4156c62ddap-1,  0x1.d906bcf328d46p-1,
+    0x1.c38b2f180bdb1p-1,  0x1.a9b66290ea1a3p-1,  0x1.8bc806b151741p-1,
+    0x1.6a09e667f3bcdp-1,  0x1.44cf325091dd6p-1,  0x1.1c73b39ae68c8p-1,
+    0x1.e2b5d3806f63bp-2,  0x1.87de2a6aea963p-2,  0x1.294062ed59f06p-2,
+    0x1.8f8b83c69a60bp-3,  0x1.917a6bc29b42cp-4,  0x0.0000000000000p+0,
+    -0x1.917a6bc29b42cp-4, -0x1.8f8b83c69a60bp-3, -0x1.294062ed59f06p-2,
+    -0x1.87de2a6aea963p-2, -0x1.e2b5d3806f63bp-2, -0x1.1c73b39ae68c8p-1,
+    -0x1.44cf325091dd6p-1, -0x1.6a09e667f3bcdp-1, -0x1.8bc806b151741p-1,
+    -0x1.a9b66290ea1a3p-1, -0x1.c38b2f180bdb1p-1, -0x1.d906bcf328d46p-1,
+    -0x1.e9f4156c62ddap-1, -0x1.f6297cff75cbp-1,  -0x1.fd88da3d12526p-1,
+    -0x1.0000000000000p+0, -0x1.fd88da3d12526p-1, -0x1.f6297cff75cbp-1,
+    -0x1.e9f4156c62ddap-1, -0x1.d906bcf328d46p-1, -0x1.c38b2f180bdb1p-1,
+    -0x1.a9b66290ea1a3p-1, -0x1.8bc806b151741p-1, -0x1.6a09e667f3bcdp-1,
+    -0x1.44cf325091dd6p-1, -0x1.1c73b39ae68c8p-1, -0x1.e2b5d3806f63bp-2,
+    -0x1.87de2a6aea963p-2, -0x1.294062ed59f06p-2, -0x1.8f8b83c69a60bp-3,
+    -0x1.917a6bc29b42cp-4,
+};
 
 static inline void sincosf_eval(double xd, uint32_t x_abs, double &sin_k,
                                 double &cos_k, double &sin_y, double &cosm1_y) {
   int64_t k;
   double y;
 
   if (likely(x_abs < FAST_PASS_BOUND)) {
     k = small_range_reduction(xd, y);
   } else {
     fputil::FPBits<float> x_bits(x_abs);
     k = large_range_reduction(xd, x_bits.get_exponent(), y);
   }
 
-  // After range reduction, k = round(x * 16 / pi) and y = (x * 16 / pi) - k.
+  // After range reduction, k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
   // So k is an integer and -0.5 <= y <= 0.5.
-  // Then sin(x) = sin((k + y)*pi/16)
-  //             = sin(y*pi/16) * cos(k*pi/16) + cos(y*pi/16) * sin(k*pi/16)
+  // Then sin(x) = sin((k + y)*pi/32)
+  //             = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
 
-  sin_k = SIN_K_PI_OVER_16[k & 31];
-  // cos(k * pi/16) = sin(k * pi/16 + pi/2) = sin((k + 8) * pi/16).
-  // cos_k = y * cos(k * pi/16)
-  cos_k = SIN_K_PI_OVER_16[(k + 8) & 31];
+  sin_k = SIN_K_PI_OVER_32[k & 63];
+  // cos(k * pi/32) = sin(k * pi/32 + pi/2) = sin((k + 16) * pi/32).
+  // cos_k = cos(k * pi/32)
+  cos_k = SIN_K_PI_OVER_32[(k + 16) & 63];
 
   double ysq = y * y;
 
-  // Degree-6 minimax even polynomial for sin(y*pi/16)/y generated by Sollya
+  // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya
   // with:
-  // > Q = fpminimax(sin(y*pi/16)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
-  sin_y = y * fputil::polyeval(ysq, 0x1.921fb54442d17p-3, -0x1.4abbce6256adp-10,
-                               0x1.466bc5a5ac6b3p-19, -0x1.32bdcb4207562p-29);
-  // Degree-8 minimax even polynomial for cos(y*pi/16) generated by Sollya with:
-  // > P = fpminimax(cos(x*pi/16), [|0, 2, 4, 6, 8|], [|1, D...|], [0, 0.5]);
-  // Note that cosm1_y = cos(y*pi/16) - 1.
-  cosm1_y =
-      ysq * fputil::polyeval(ysq, -0x1.3bd3cc9be45dcp-6, 0x1.03c1f081b08ap-14,
-                             -0x1.55d3c6fb0fb6ep-24, 0x1.e1d3d60f58873p-35);
+  // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
+  sin_y =
+      y * fputil::polyeval(ysq, 0x1.921fb54442d18p-4, -0x1.4abbce625abb1p-13,
+                           0x1.466bc624f2776p-24, -0x1.32c3a619d4a7ep-36);
+  // Degree-8 minimax even polynomial for cos(y*pi/32) generated by Sollya with:
+  // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, D...|], [0, 0.5]);
+  // Note that cosm1_y = cos(y*pi/32) - 1.
+  cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3cc9be430bp-8,
+                                   0x1.03c1f070c2e27p-18, -0x1.55cc84bd942p-30);
 }
 
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H

libc/src/math/generic/sinf.cpp

 //===-- Single-precision sin function -------------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "src/math/sinf.h"
 #include "sincosf_utils.h"
 #include "src/__support/FPUtil/BasicOperations.h"
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PolyEval.h"
-#include "src/__support/FPUtil/except_value_utils.h"
 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/common.h"
 
 #include <errno.h>
 
 #if defined(LIBC_TARGET_HAS_FMA)
 #include "range_reduction_fma.h"
-// using namespace __llvm_libc::fma;
-using __llvm_libc::fma::N_EXCEPTS;
-using __llvm_libc::fma::SinfExcepts;
 #else
 #include "range_reduction.h"
-// using namespace __llvm_libc::generic;
-using __llvm_libc::generic::N_EXCEPTS;
-using __llvm_libc::generic::SinfExcepts;
 #endif
 
 namespace __llvm_libc {
 
 LLVM_LIBC_FUNCTION(float, sinf, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
 
   uint32_t x_u = xbits.uintval();
   uint32_t x_abs = x_u & 0x7fff'ffffU;
   double xd = static_cast<double>(x);
 
   // Range reduction:
-  // For |x| > pi/16, we perform range reduction as follows:
+  // For |x| > pi/32, we perform range reduction as follows:
   // Find k and y such that:
-  //   x = (k + y) * pi/16
+  //   x = (k + y) * pi/32
   //   k is an integer
   //   |y| < 0.5
-  // For small range (|x| < 2^46 when FMA instructions are available, 2^22
+  // For small range (|x| < 2^45 when FMA instructions are available, 2^22
   // otherwise), this is done by performing:
-  //   k = round(x * 16/pi)
-  //   y = x * 16/pi - k
-  // For large range, we will omit all the higher parts of 16/pi such that the
-  // least significant bits of their full products with x are larger than 31,
-  // since sin((k + y + 32*i) * pi/16) = sin(x + i * 2pi) = sin(x).
+  //   k = round(x * 32/pi)
+  //   y = x * 32/pi - k
+  // For large range, we will omit all the higher parts of 32/pi such that the
+  // least significant bits of their full products with x are larger than 63,
+  // since sin((k + y + 64*i) * pi/32) = sin(x + i * 2pi) = sin(x).
   //
-  // When FMA instructions are not available, we store the digits of 16/pi in
+  // When FMA instructions are not available, we store the digits of 32/pi in
   // chunks of 28-bit precision.  This will make sure that the products:
-  //   x * SIXTEEN_OVER_PI_28[i] are all exact.
-  // When FMA instructions are available, we simply store the digits of 16/pi in
+  //   x * THIRTYTWO_OVER_PI_28[i] are all exact.
+  // When FMA instructions are available, we simply store the digits of 32/pi in
   // chunks of doubles (53-bit of precision).
   // So when multiplying by the largest values of single precision, the
   // resulting output should be correct up to 2^(-208 + 128) ~ 2^-80.  By the
   // worst-case analysis of range reduction, |y| >= 2^-38, so this should give
   // us more than 40 bits of accuracy. For the worst-case estimation of range
   // reduction, see for instances:
   //   Elementary Functions by J-M. Muller, Chapter 11,
   //   Handbook of Floating-Point Arithmetic by J-M. Muller et. al.,
   //   Chapter 10.2.
   //
   // Once k and y are computed, we then deduce the answer by the sine of sum
   // formula:
-  //   sin(x) = sin((k + y)*pi/16)
-  //          = sin(y*pi/16) * cos(k*pi/16) + cos(y*pi/16) * sin(k*pi/16)
-  // The values of sin(k*pi/16) and cos(k*pi/16) for k = 0..31 are precomputed
-  // and stored using a vector of 32 doubles. Sin(y*pi/16) and cos(y*pi/16) are
-  // computed using degree-7 and degree-8 minimax polynomials generated by
+  //   sin(x) = sin((k + y)*pi/32)
+  //          = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
+  // The values of sin(k*pi/32) and cos(k*pi/32) for k = 0..31 are precomputed
+  // and stored using a vector of 32 doubles. Sin(y*pi/32) and cos(y*pi/32) are
+  // computed using degree-7 and degree-6 minimax polynomials generated by
   // Sollya respectively.
 
   // |x| <= pi/16
   if (unlikely(x_abs <= 0x3e49'0fdbU)) {
 
     // |x| < 0x1.d12ed2p-12f
     if (unlikely(x_abs < 0x39e8'9769U)) {
       if (unlikely(x_abs == 0U)) {
#endif // LIBC_TARGET_HAS_FMA
     // > display = hexadecimal;
     // > Q = fpminimax(sin(x)/x, [|0, 2, 4, 6, 8|], [|1, D...|], [0, pi/16]);
     double result =
         fputil::polyeval(xsq, 1.0, -0x1.55555555554c6p-3, 0x1.1111111085e65p-7,
                          -0x1.a019f70fb4d4fp-13, 0x1.718d179815e74p-19);
     return xd * result;
   }
 
-  using ExceptChecker = typename fputil::ExceptionChecker<float, N_EXCEPTS>;
-  {
-    float result;
-    if (ExceptChecker::check_odd_func(SinfExcepts, x_abs, xbits.get_sign(),
-                                      result))
-      return result;
+  if (unlikely(x_abs == 0x4619'9998U)) { // x = 0x1.33333p13
+    float r = -0x1.63f4bap-2f;
+    int rounding = fputil::get_round();
+    bool sign = xbits.get_sign();
+    if ((rounding == FE_DOWNWARD && !sign) || (rounding == FE_UPWARD && sign))
+      r = -0x1.63f4bcp-2f;
+    return xbits.get_sign() ? -r : r;
   }
 
   if (unlikely(x_abs >= 0x7f80'0000U)) {
     if (x_abs == 0x7f80'0000U) {
       errno = EDOM;
       fputil::set_except(FE_INVALID);
     }
     return x +
            FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
   }
 
   // Combine the results with the sine of sum formula:
-  //   sin(x) = sin((k + y)*pi/16)
-  //          = sin(y*pi/16) * cos(k*pi/16) + cos(y*pi/16) * sin(k*pi/16)
+  //   sin(x) = sin((k + y)*pi/32)
+  //          = sin(y*pi/32) * cos(k*pi/32) + cos(y*pi/32) * sin(k*pi/32)
   //          = sin_y * cos_k + (1 + cosm1_y) * sin_k
   //          = sin_y * cos_k + (cosm1_y * sin_k + sin_k)
   double sin_k, cos_k, sin_y, cosm1_y;
 
   sincosf_eval(xd, x_abs, sin_k, cos_k, sin_y, cosm1_y);
 
   return fputil::multiply_add(sin_y, cos_k,
                               fputil::multiply_add(cosm1_y, sin_k, sin_k));
 }
 
 } // namespace __llvm_libc