[libc] Implement sincosf function correctly rounded to all rounding modes.

Refactor common range reductions and evaluations for sinf, cosf, and
sincosf.  Added exhaustive tests for sincosf.

Performance before the patch:
```
System LIBC reciprocal throughput : 30.205
LIBC reciprocal throughput        : 30.533

System LIBC latency : 67.961
LIBC latency        : 61.564
```
Performance after the patch:
```
System LIBC reciprocal throughput : 30.409
LIBC reciprocal throughput        : 20.273

System LIBC latency : 67.527
LIBC latency        : 61.959
```

Reviewed By: orex

Differential Revision: https://reviews.llvm.org/D130901

Author: lntue

libc/docs/math.rst

@@ -166,7 +166,7 @@ log10          |check|
 log1p          |check|
 log2           |check|
 sin            |check|          large
-sincos         0.776 ULPs       large
+sincos         |check|          large
 sinh           |check|
 sqrt           |check|          |check|         |check|
 tanh           |check|
@@ -232,6 +232,8 @@ Performance
 +--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | sinf         |        13 |                25 |        54 |                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           |
++--------------+-----------+-------------------+-----------+-------------------+-------------------------------------+------------+-------------------------+--------------+---------------+
 | 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           |

libc/src/math/generic/CMakeLists.txt

@@ -46,14 +46,26 @@ add_object_library(
     libc.src.errno.errno
 )
 
-add_object_library(
+add_header_library(
+  range_reduction
+  HDRS
+    range_reduction.h
+    range_reduction_fma.h
+  DEPENDS
+    libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.fma
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.nearest_integer
+)
+
+add_header_library(
   sincosf_utils
   HDRS
     sincosf_utils.h
-  SRCS
-    sincosf_data.cpp
   DEPENDS
-    .math_utils
+    .range_reduction
+    libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.polyeval
 )
 
 add_entrypoint_object(
@@ -62,16 +74,13 @@ add_entrypoint_object(
     cosf.cpp
   HDRS
     ../cosf.h
-    range_reduction.h
-    range_reduction_fma.h
   DEPENDS
-    .common_constants
+    .sincosf_utils
     libc.include.math
     libc.src.errno.errno
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.fma
     libc.src.__support.FPUtil.multiply_add
-    libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
   COMPILE_OPTIONS
     -O3
@@ -83,16 +92,14 @@ add_entrypoint_object(
     sinf.cpp
   HDRS
     ../sinf.h
-    range_reduction.h
-    range_reduction_fma.h
   DEPENDS
-    .common_constants
+    .range_reduction
+    .sincosf_utils
     libc.include.math
     libc.src.errno.errno
     libc.src.__support.FPUtil.fputil
     libc.src.__support.FPUtil.fma
     libc.src.__support.FPUtil.multiply_add
-    libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
   COMPILE_OPTIONS
     -O3
@@ -105,9 +112,14 @@ add_entrypoint_object(
   HDRS
     ../sincosf.h
   DEPENDS
+    .range_reduction
     .sincosf_utils
     libc.include.math
     libc.src.errno.errno
+    libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.fma
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.polyeval
   COMPILE_OPTIONS
     -O3
 )

libc/src/math/generic/common_constants.cpp

@@ -240,21 +240,4 @@ const double EXP_2_POW[EXP_num_p] = {
     0x1.13821818624b4p-2,  0x1.2ff6b54d8a89cp-2,  0x1.4d0ad5a753e07p-2,
     0x1.6ac1f752150a5p-2,  0x1.891fac0e95613p-2};
 
-// Lookup table for sin(k * pi / 16) with k = 0, ..., 31.
-// 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};
-
 } // namespace __llvm_libc

libc/src/math/generic/common_constants.h

@@ -31,18 +31,13 @@ extern const double EXP_M1[195];
 // > for i from 0 to 127 do { D(exp(i / 128)); };
 extern const double EXP_M2[128];
 
-// Lookup table for sin(k * pi / 16) with k = 0, ..., 31.
-// Table is generated with Sollya as follow:
-// > display = hexadecimal;
-// > for k from 0 to 31 do { D(sin(k * pi/16)); };
-extern const double SIN_K_PI_OVER_16[32];
-
 static constexpr int EXP_bits_p = 5;
 static constexpr int EXP_num_p = 1 << EXP_bits_p;
 
 // Wolfram alpha: N[Table[2^x-1,{x,-16/32,15/32,1/32}],27]
 // printf("%.13a,\n", d[i]);
 extern const double EXP_2_POW[EXP_num_p];
+
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_MATH_GENERIC_COMMON_CONSTANTS_H

libc/src/math/generic/cosf.cpp

@@ -7,31 +7,16 @@
 //===----------------------------------------------------------------------===//
 
 #include "src/math/cosf.h"
-#include "common_constants.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::FAST_PASS_BOUND;
-using __llvm_libc::fma::large_range_reduction;
-using __llvm_libc::fma::small_range_reduction;
-#else
-#include "range_reduction.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
-
 namespace __llvm_libc {
 
 // Exceptional cases for cosf.
@@ -132,58 +117,26 @@ LLVM_LIBC_FUNCTION(float, cosf, (float x)) {
       return result;
   }
 
-  // TODO(lntue): refactor range reduction and most of polynomial approximation
-  // to share between sinf, cosf, and sincosf.
-  int k;
-  double y;
-
-  if (likely(x_abs < FAST_PASS_BOUND)) {
-    k = small_range_reduction(xd, y);
-  } else {
-    // x is inf or nan.
-    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));
+  // x is inf or nan.
+  if (unlikely(x_abs >= 0x7f80'0000U)) {
+    if (x_abs == 0x7f80'0000U) {
+      errno = EDOM;
+      fputil::set_except(FE_INVALID);
     }
-
-    k = large_range_reduction(xd, xbits.get_exponent(), y);
+    return x +
+           FPBits::build_nan(1 << (fputil::MantissaWidth<float>::VALUE - 1));
   }
 
-  // After range reduction, k = round(x * 16 / pi) and y = (x * 16 / pi) - k.
-  // So k is an integer and -0.5 <= y <= 0.5.
-  // Then cos(x) = cos((k + y)*pi/16)
-  //             = cos(y*pi/16) * cos(k*pi/16) - sin(y*pi/16) * sin(k*pi/16)
-
-  double ysq = y * y;
-
-  // Degree-6 minimax even polynomial for sin(y*pi/16)/y generated by Sollya
-  // with:
-  // > Q = fpminimax(sin(y*pi/16)/y, [|0, 2, 4, 6|], [|D...|], [0, 0.5]);
-  double 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.
-  double cosm1_y =
-      ysq * fputil::polyeval(ysq, -0x1.3bd3cc9be45dcp-6, 0x1.03c1f081b08ap-14,
-                             -0x1.55d3c6fb0fb6ep-24, 0x1.e1d3d60f58873p-35);
-
-  double 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)
-  double cos_k = SIN_K_PI_OVER_16[(k + 8) & 31];
-
   // 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)
   //          = cosm1_y * cos_k + sin_y * sin_k
   //          = (cosm1_y * cos_k + cos_k) + sin_y * sin_k
-  return fputil::multiply_add(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));
 }