[libc][math] Improve performance of double precision trig functions. (#111793)

- Improve the accuracy of fast pass' range reduction.
- Provide tighter error estimations.
- Reduce the table size when `LIBC_MATH_SMALL_TABLES` flag is set.

Author: lntue

libc/src/__support/FPUtil/double_double.h

@@ -18,6 +18,8 @@
 namespace LIBC_NAMESPACE_DECL {
 namespace fputil {
 
+#define DEFAULT_DOUBLE_SPLIT 27
+
 using DoubleDouble = LIBC_NAMESPACE::NumberPair<double>;
 
 // The output of Dekker's FastTwoSum algorithm is correct, i.e.:
@@ -61,7 +63,8 @@ LIBC_INLINE constexpr DoubleDouble add(const DoubleDouble &a, double b) {
 //   Zimmermann, P., "Note on the Veltkamp/Dekker Algorithms with Directed
 //   Roundings," https://inria.hal.science/hal-04480440.
 // Default splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1.
-template <size_t N = 27> LIBC_INLINE constexpr DoubleDouble split(double a) {
+template <size_t N = DEFAULT_DOUBLE_SPLIT>
+LIBC_INLINE constexpr DoubleDouble split(double a) {
   DoubleDouble r{0.0, 0.0};
   // CN = 2^N.
   constexpr double CN = static_cast<double>(1 << N);
@@ -73,14 +76,30 @@ template <size_t N = 27> LIBC_INLINE constexpr DoubleDouble split(double a) {
   return r;
 }
 
+// Helper for non-fma exact mult where the first number is already split.
+template <size_t SPLIT_B = DEFAULT_DOUBLE_SPLIT>
+LIBC_INLINE DoubleDouble exact_mult(const DoubleDouble &as, double a,
+                                    double b) {
+  DoubleDouble bs = split<SPLIT_B>(b);
+  DoubleDouble r{0.0, 0.0};
+
+  r.hi = a * b;
+  double t1 = as.hi * bs.hi - r.hi;
+  double t2 = as.hi * bs.lo + t1;
+  double t3 = as.lo * bs.hi + t2;
+  r.lo = as.lo * bs.lo + t3;
+
+  return r;
+}
+
 // Note: When FMA instruction is not available, the `exact_mult` function is
 // only correct for round-to-nearest mode.  See:
 //   Zimmermann, P., "Note on the Veltkamp/Dekker Algorithms with Directed
 //   Roundings," https://inria.hal.science/hal-04480440.
 // Using Theorem 1 in the paper above, without FMA instruction, if we restrict
 // the generated constants to precision <= 51, and splitting it by 2^28 + 1,
 // then a * b = r.hi + r.lo is exact for all rounding modes.
-template <bool NO_FMA_ALL_ROUNDINGS = false>
+template <size_t SPLIT_B = 27>
 LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
   DoubleDouble r{0.0, 0.0};
 
@@ -90,18 +109,8 @@ LIBC_INLINE DoubleDouble exact_mult(double a, double b) {
 #else
   // Dekker's Product.
   DoubleDouble as = split(a);
-  DoubleDouble bs;
 
-  if constexpr (NO_FMA_ALL_ROUNDINGS)
-    bs = split<28>(b);
-  else
-    bs = split(b);
-
-  r.hi = a * b;
-  double t1 = as.hi * bs.hi - r.hi;
-  double t2 = as.hi * bs.lo + t1;
-  double t3 = as.lo * bs.hi + t2;
-  r.lo = as.lo * bs.lo + t3;
+  r = exact_mult<SPLIT_B>(as, a, b);
 #endif // LIBC_TARGET_CPU_HAS_FMA
 
   return r;
@@ -113,10 +122,10 @@ LIBC_INLINE DoubleDouble quick_mult(double a, const DoubleDouble &b) {
   return r;
 }
 
-template <bool NO_FMA_ALL_ROUNDINGS = false>
+template <size_t SPLIT_B = 27>
 LIBC_INLINE DoubleDouble quick_mult(const DoubleDouble &a,
                                     const DoubleDouble &b) {
-  DoubleDouble r = exact_mult<NO_FMA_ALL_ROUNDINGS>(a.hi, b.hi);
+  DoubleDouble r = exact_mult<SPLIT_B>(a.hi, b.hi);
   double t1 = multiply_add(a.hi, b.lo, r.lo);
   double t2 = multiply_add(a.lo, b.hi, t1);
   r.lo = t2;
@@ -157,8 +166,8 @@ LIBC_INLINE DoubleDouble div(const DoubleDouble &a, const DoubleDouble &b) {
   double e_hi = fputil::multiply_add(b.hi, -r.hi, a.hi);
   double e_lo = fputil::multiply_add(b.lo, -r.hi, a.lo);
 #else
-  DoubleDouble b_hi_r_hi = fputil::exact_mult</*NO_FMA=*/true>(b.hi, -r.hi);
-  DoubleDouble b_lo_r_hi = fputil::exact_mult</*NO_FMA=*/true>(b.lo, -r.hi);
+  DoubleDouble b_hi_r_hi = fputil::exact_mult(b.hi, -r.hi);
+  DoubleDouble b_lo_r_hi = fputil::exact_mult(b.lo, -r.hi);
   double e_hi = (a.hi + b_hi_r_hi.hi) + b_hi_r_hi.lo;
   double e_lo = (a.lo + b_lo_r_hi.hi) + b_lo_r_hi.lo;
 #endif // LIBC_TARGET_CPU_HAS_FMA

libc/src/__support/macros/optimization.h

@@ -48,6 +48,16 @@ LIBC_INLINE constexpr bool expects_bool_condition(T value, T expected) {
 
 #ifndef LIBC_MATH
 #define LIBC_MATH 0
+#else
+
+#if (LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS)
+#define LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#endif
+
+#if (LIBC_MATH & LIBC_MATH_SMALL_TABLES)
+#define LIBC_MATH_HAS_SMALL_TABLES
+#endif
+
 #endif // LIBC_MATH
 
 #endif // LLVM_LIBC_SRC___SUPPORT_MACROS_OPTIMIZATION_H

libc/src/math/generic/cos.cpp

@@ -17,17 +17,14 @@
 #include "src/__support/macros/config.h"
 #include "src/__support/macros/optimization.h"            // LIBC_UNLIKELY
 #include "src/__support/macros/properties/cpu_features.h" // LIBC_TARGET_CPU_HAS_FMA
+#include "src/math/generic/range_reduction_double_common.h"
 #include "src/math/generic/sincos_eval.h"
 
-// TODO: We might be able to improve the performance of large range reduction of
-// non-FMA targets further by operating directly on 25-bit chunks of 128/pi and
-// pre-split SIN_K_PI_OVER_128, but that might double the memory footprint of
-// those lookup table.
-#include "range_reduction_double_common.h"
-
-#if ((LIBC_MATH & LIBC_MATH_SKIP_ACCURATE_PASS) != 0)
-#define LIBC_MATH_COS_SKIP_ACCURATE_PASS
-#endif
+#ifdef LIBC_TARGET_CPU_HAS_FMA
+#include "range_reduction_double_fma.h"
+#else
+#include "range_reduction_double_nofma.h"
+#endif // LIBC_TARGET_CPU_HAS_FMA
 
 namespace LIBC_NAMESPACE_DECL {
 
@@ -42,22 +39,29 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) {
 
   DoubleDouble y;
   unsigned k;
-  generic::LargeRangeReduction<NO_FMA> range_reduction_large{};
+  LargeRangeReduction range_reduction_large{};
 
-  // |x| < 2^32 (with FMA) or |x| < 2^23 (w/o FMA)
+  // |x| < 2^16.
   if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT)) {
-    // |x| < 2^-27
-    if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) {
-      // Signed zeros.
-      if (LIBC_UNLIKELY(x == 0.0))
-        return 1.0;
-
-      // For |x| < 2^-27, |cos(x) - 1| < |x|^2/2 < 2^-54 = ulp(1 - 2^-53)/2.
-      return fputil::round_result_slightly_down(1.0);
+    // |x| < 2^-7
+    if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 7)) {
+      // |x| < 2^-27
+      if (LIBC_UNLIKELY(x_e < FPBits::EXP_BIAS - 27)) {
+        // Signed zeros.
+        if (LIBC_UNLIKELY(x == 0.0))
+          return 1.0;
+
+        // For |x| < 2^-27, |cos(x) - 1| < |x|^2/2 < 2^-54 = ulp(1 - 2^-53)/2.
+        return fputil::round_result_slightly_down(1.0);
+      }
+      // No range reduction needed.
+      k = 0;
+      y.lo = 0.0;
+      y.hi = x;
+    } else {
+      // Small range reduction.
+      k = range_reduction_small(x, y);
     }
-
-    // // Small range reduction.
-    k = range_reduction_small(x, y);
   } else {
     // Inf or NaN
     if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
@@ -70,70 +74,51 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) {
     }
 
     // Large range reduction.
-    k = range_reduction_large.compute_high_part(x);
-    y = range_reduction_large.fast();
+    k = range_reduction_large.fast(x, y);
   }
 
   DoubleDouble sin_y, cos_y;
 
-  generic::sincos_eval(y, sin_y, cos_y);
+  [[maybe_unused]] double err = generic::sincos_eval(y, sin_y, cos_y);
 
   // Look up sin(k * pi/128) and cos(k * pi/128)
-  // Memory saving versions:
-
-  // Use 128-entry table instead:
-  // DoubleDouble sin_k = SIN_K_PI_OVER_128[k & 127];
-  // uint64_t sin_s = static_cast<uint64_t>((k + 128) & 128) << (63 - 7);
-  // sin_k.hi = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val();
-  // sin_k.lo = FPBits(FPBits(sin_k.hi).uintval() ^ sin_s).get_val();
-  // DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 127];
-  // uint64_t cos_s = static_cast<uint64_t>((k + 64) & 128) << (63 - 7);
-  // cos_k.hi = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val();
-  // cos_k.lo = FPBits(FPBits(cos_k.hi).uintval() ^ cos_s).get_val();
-
-  // Use 64-entry table instead:
-  // auto get_idx_dd = [](unsigned kk) -> DoubleDouble {
-  //   unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
-  //   DoubleDouble ans = SIN_K_PI_OVER_128[idx];
-  //   if (kk & 128) {
-  //     ans.hi = -ans.hi;
-  //     ans.lo = -ans.lo;
-  //   }
-  //   return ans;
-  // };
-  // DoubleDouble sin_k = get_idx_dd(k + 128);
-  // DoubleDouble cos_k = get_idx_dd(k + 64);
-
+#ifdef LIBC_MATH_HAS_SMALL_TABLES
+  // Memory saving versions.  Use 65-entry table.
+  auto get_idx_dd = [](unsigned kk) -> DoubleDouble {
+    unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
+    DoubleDouble ans = SIN_K_PI_OVER_128[idx];
+    if (kk & 128) {
+      ans.hi = -ans.hi;
+      ans.lo = -ans.lo;
+    }
+    return ans;
+  };
+  DoubleDouble sin_k = get_idx_dd(k + 128);
+  DoubleDouble cos_k = get_idx_dd(k + 64);
+#else
   // Fast look up version, but needs 256-entry table.
   // -sin(k * pi/128) = sin((k + 128) * pi/128)
   // cos(k * pi/128) = sin(k * pi/128 + pi/2) = sin((k + 64) * pi/128).
   DoubleDouble msin_k = SIN_K_PI_OVER_128[(k + 128) & 255];
   DoubleDouble cos_k = SIN_K_PI_OVER_128[(k + 64) & 255];
+#endif // LIBC_MATH_HAS_SMALL_TABLES
 
   // After range reduction, k = round(x * 128 / pi) and y = x - k * (pi / 128).
   // So k is an integer and -pi / 256 <= y <= pi / 256.
   // Then cos(x) = cos((k * pi/128 + y)
   //             = cos(y) * cos(k*pi/128) - sin(y) * sin(k*pi/128)
-  DoubleDouble cos_k_cos_y = fputil::quick_mult<NO_FMA>(cos_y, cos_k);
-  DoubleDouble msin_k_sin_y = fputil::quick_mult<NO_FMA>(sin_y, msin_k);
+  DoubleDouble cos_k_cos_y = fputil::quick_mult(cos_y, cos_k);
+  DoubleDouble msin_k_sin_y = fputil::quick_mult(sin_y, msin_k);
 
   DoubleDouble rr = fputil::exact_add<false>(cos_k_cos_y.hi, msin_k_sin_y.hi);
   rr.lo += msin_k_sin_y.lo + cos_k_cos_y.lo;
 
-#ifdef LIBC_MATH_COS_SKIP_ACCURATE_PASS
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
   return rr.hi + rr.lo;
 #else
 
-  // Accurate test and pass for correctly rounded implementation.
-#ifdef LIBC_TARGET_CPU_HAS_FMA
-  constexpr double ERR = 0x1.0p-70;
-#else
-  // TODO: Improve non-FMA fast pass accuracy.
-  constexpr double ERR = 0x1.0p-66;
-#endif // LIBC_TARGET_CPU_HAS_FMA
-
-  double rlp = rr.lo + ERR;
-  double rlm = rr.lo - ERR;
+  double rlp = rr.lo + err;
+  double rlm = rr.lo - err;
 
   double r_upper = rr.hi + rlp; // (rr.lo + ERR);
   double r_lower = rr.hi + rlm; // (rr.lo - ERR);
@@ -144,15 +129,15 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) {
 
   Float128 u_f128, sin_u, cos_u;
   if (LIBC_LIKELY(x_e < FPBits::EXP_BIAS + FAST_PASS_EXPONENT))
-    u_f128 = generic::range_reduction_small_f128(x);
+    u_f128 = range_reduction_small_f128(x);
   else
     u_f128 = range_reduction_large.accurate();
 
   generic::sincos_eval(u_f128, sin_u, cos_u);
 
   auto get_sin_k = [](unsigned kk) -> Float128 {
     unsigned idx = (kk & 64) ? 64 - (kk & 63) : (kk & 63);
-    Float128 ans = generic::SIN_K_PI_OVER_128_F128[idx];
+    Float128 ans = SIN_K_PI_OVER_128_F128[idx];
     if (kk & 128)
       ans.sign = Sign::NEG;
     return ans;
@@ -172,7 +157,7 @@ LLVM_LIBC_FUNCTION(double, cos, (double x)) {
   // https://github.com/llvm/llvm-project/issues/96452.
 
   return static_cast<double>(r);
-#endif // !LIBC_MATH_COS_SKIP_ACCURATE_PASS
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
 }
 
 } // namespace LIBC_NAMESPACE_DECL

libc/src/math/generic/pow.cpp

@@ -398,7 +398,7 @@ LLVM_LIBC_FUNCTION(double, pow, (double x, double y)) {
 #else
   double c = FPBits(m_x.uintval() & 0x3fff'e000'0000'0000).get_val();
   dx = fputil::multiply_add(RD[idx_x], m_x.get_val() - c, CD[idx_x]); // Exact
-  dx_c0 = fputil::exact_mult<true>(COEFFS[0], dx);
+  dx_c0 = fputil::exact_mult<28>(dx, COEFFS[0]);                      // Exact
 #endif // LIBC_TARGET_CPU_HAS_FMA
 
   double dx2 = dx * dx;