BUG: Fix f64 trig accuracy; match SVML rounding DAGs

npsr/hwy.h

@@ -1,7 +1,7 @@
 #ifndef NPSR_HWY_H_
 #define NPSR_HWY_H_
 
 #include <hwy/highway.h>
 // This macro is used to define intrinsics that are:
 // NOTE: equals to HWY_API.
 // - always inlined
@@ -21,7 +21,7 @@
 #endif
 
 HWY_BEFORE_NAMESPACE();
 namespace npsr::HWY_NAMESPACE {
 namespace hn = hwy::HWY_NAMESPACE;
 using hn::DFromV;
 using hn::MFromD;
@@ -32,7 +32,23 @@
 using hn::VFromD;
 constexpr bool kNativeFMA = HWY_NATIVE_FMA != 0;
 
-inline HWY_ATTR void DummyToSuppressUnusedWarning() {}
+// Dekker split product for targets without native FMA, where the
+// `MulSub(a, b, a*b)` error idiom degenerates to zero. Masking to 26-bit heads
+// makes `head` exact, so a*b == head + rest (rest rounded at ~2^-105).
+template <typename V, HWY_IF_F64(TFromV<V>)>
+NPSR_INTRIN void SplitMul(V a, V b, V& head, V& rest) {
+  using namespace hn;
+  const DFromV<V> d;
+  const RebindToUnsigned<decltype(d)> du;
+  const V mask = BitCast(d, Set(du, 0xFFFFFFFFF8000000u));
+  const V a_hi = And(a, mask);
+  const V a_lo = Sub(a, a_hi);
+  const V b_hi = And(b, mask);
+  const V b_lo = Sub(b, b_hi);
+  head = Mul(a_hi, b_hi);
+  rest = MulAdd(a_hi, b_lo, Mul(a_lo, b));
+}
+
 }  // namespace npsr::HWY_NAMESPACE
 HWY_AFTER_NAMESPACE();
 

npsr/lut-inl.h

@@ -1,3 +1,3 @@
 #if defined(NPSR_LUT_INL_H_) == defined(HWY_TARGET_TOGGLE)  // NOLINT
 #ifdef NPSR_LUT_INL_H_
 #undef NPSR_LUT_INL_H_
@@ -5,13 +5,14 @@
 #define NPSR_LUT_INL_H_
 #endif
 
+#include <limits>
 #include <tuple>
 
 #include "npsr/hwy.h"
 
 HWY_BEFORE_NAMESPACE();
 
 namespace npsr::HWY_NAMESPACE {
 
 /**
  * @brief Optimized Lookup Table.

npsr/trig/data/constants.h

@@ -1,9 +1,9 @@
 // Auto-generated by npsr/trig/data/constants.h.sol
 // Use `spin sollya -f` to force regeneration
 #ifndef NPSR_TRIG_DATA_CONSTANTS_H
 #define NPSR_TRIG_DATA_CONSTANTS_H
 
 namespace npsr::trig::data {
 template <typename T, bool FMA> inline constexpr char kPi[] = {};
 template <> inline constexpr float kPi<float, true>[] = {
 0x1.921fb6p1f, -0x1.777a5cp-24f, -0x1.ee59dap-49f, };
@@ -30,7 +30,7 @@
 template <bool FMA> inline constexpr double kPiDiv16Prec29[] = {
 0x1.921fb54442d18p-3, 0x1.1a62633p-57, 0x1.45c06e0e68948p-89, };
 template <> inline constexpr double kPiDiv16Prec29<false>[] = {
-0x1.921fb54p-3, 0x1.a626331p-61, 0x1.1701b839a252p-91, 0x1.10b461p-33, 
+0x1.921fb54p-3, 0x1.10b461p-33, 0x1.a626331p-61, 0x1.1701b839a252p-91, 
 };
 
 template <typename T> inline constexpr char kInvPi = '_';

npsr/trig/extended-inl.h

@@ -5,12 +5,12 @@
 #define NPSR_TRIG_EXTENDED_INL_H_
 #endif
 
 #include "npsr/hwy.h"
 #include "npsr/trig/data/data.h"
 #include "npsr/trig/low-inl.h"  // Operation
 
 HWY_BEFORE_NAMESPACE();
 namespace npsr::HWY_NAMESPACE::trig {
 
 template <Operation OP, class V>
 NPSR_INTRIN V Extended(V x) {
@@ -212,7 +212,7 @@
   // PHASE 7: Convert to Radians
   // =============================================================================
 
-  // Multiply by π with error compensation (Cody-Waite multiplication)
+  // Multiply by 2π with error compensation (Cody-Waite multiplication)
   constexpr auto kPiMul2 = data::kPiMul2<T>;
   const V pi2_hi = Set(d, kPiMul2[0]);
   const V pi2_med = Set(d, kPiMul2[1]);
@@ -233,6 +233,14 @@
     r_lo = Combine(d, DemoteTo(dh, r_lo_w0), DemoteTo(dh, r_lo_w1));
     r_w0 = BitCast(d, r0);
     r_w1 = BitCast(d, r1);
+  } else if constexpr (!kNativeFMA) {
+    // Without FMA, `MulSub(pi2_hi, n, r)` yields zero and drops the 0.5-ulp(r)
+    // rounding error of r — worth 0.5 ulp of the result near k*pi/2. Recover
+    // it from the exact head product (head - r is exact by Sterbenz).
+    V head, rest;
+    SplitMul(pi2_hi, n, head, rest);
+    r_lo = Add(Sub(head, r), rest);
+    r_lo = Add(Mul(pi2_med, n), r_lo);
   } else {
     r_lo = MulSub(pi2_hi, n, r);
     r_lo = MulAdd(pi2_med, n, r_lo);
@@ -272,22 +280,30 @@
   // =============================================================================
   // PHASE 10: Final Assembly
   // =============================================================================
-  V res_lo = NegMulAdd(func_hi, r, deriv);
-  res_lo = MulAdd(res_lo, r_lo, func_lo);
   V res_hi_lo = MulAdd(sigma, r, func_hi);
-  V res_hi = MulAdd(deriv_hi, r, res_hi_lo);
-
+  V res_hi, deriv_hi_r_cor;
+  if constexpr (!kNativeFMA && !kIsSingle) {
+    // Same substitution as in Phase 7. sigma is a signed power of two (see
+    // approx.h) so sigma*r is exact, but deriv_hi*r would round under a
+    // decayed MulAdd and corrupt the head; SplitMul recovers it.
+    V dr, dr_rest;
+    SplitMul(deriv_hi, r, dr, dr_rest);
+    res_hi = Add(res_hi_lo, dr);
+    V dr_hi = Sub(res_hi, res_hi_lo);
+    deriv_hi_r_cor = Add(Sub(dr, dr_hi), dr_rest);
+  } else {
+    res_hi = MulAdd(deriv_hi, r, res_hi_lo);
+    deriv_hi_r_cor = MulAdd(deriv_hi, r, Sub(res_hi_lo, res_hi));
+  }
   V sum_cor = MulAdd(sigma, r, Sub(func_hi, res_hi_lo));
-  V deriv_hi_r_cor = MulAdd(deriv_hi, r, Sub(res_hi_lo, res_hi));
-  deriv_hi_r_cor = Add(deriv_hi_r_cor, sum_cor);
-  res_lo = Add(res_lo, deriv_hi_r_cor);
+  V res_lo0 = Add(deriv_hi_r_cor, sum_cor);
 
   // Polynomial corrections
   V s2 = Set(d, kIsSingle ? 0x1.1110b8p-7f : 0x1.1110fabb3551cp-7);
   V s1 = Set(d, kIsSingle ? -0x1.555556p-3f : -0x1.5555555554448p-3);
   V sin_poly = MulAdd(s2, r2, s1);
-  sin_poly = Mul(sin_poly, r);
   sin_poly = Mul(sin_poly, r2);
+  sin_poly = Mul(sin_poly, r);
 
   V c1 = Set(d, kIsSingle ? 0x1.5554f8p-5f : 0x1.5555555554ccfp-5);
   const V neg_half = Set(d, static_cast<T>(-0.5));
@@ -301,8 +317,11 @@
   }
   cos_poly = Mul(cos_poly, r2);
 
-  res_lo = MulAdd(sin_poly, deriv, res_lo);
-  res_lo = MulAdd(cos_poly, func_hi, res_lo);
+  V corr = NegMulAdd(func_hi, r, deriv);
+  V res_lo_rlo = MulAdd(corr, r_lo, func_lo);
+  V res_lo_sin = MulAdd(deriv, sin_poly, res_lo0);
+  V res_lo_cos = MulAdd(func_hi, cos_poly, res_lo_sin);
+  V res_lo = Add(res_lo_cos, res_lo_rlo);
   return Add(res_hi, res_lo);
 }
 // NOLINTNEXTLINE(google-readability-namespace-comments)

npsr/trig/high-inl.h

@@ -5,14 +5,14 @@
 #define NPSR_TRIG_HIGH_INL_H_
 #endif
 
 #include "npsr/hwy.h"
 #include "npsr/lut-inl.h"
 #include "npsr/trig/data/data.h"
 #include "npsr/trig/low-inl.h"  // Operation
 
 HWY_BEFORE_NAMESPACE();
 
 namespace npsr::HWY_NAMESPACE::trig {
 
 template <Operation OP, typename V, HWY_IF_F32(TFromV<V>)>
 NPSR_INTRIN V High(V x) {
@@ -95,15 +95,16 @@
 /**
  * This function computes sin(x) or cos(x) for |x| < 2^24 using the Cody-Waite
  * reduction algorithm combined with table lookup and polynomial approximation,
- * achieves < 1 ULP error for |x| < 2^24.
+ * achieves <= 1 ULP error for |x| < 2^24 (worst case reaches, but does not
+ * exceed, 1 ULP).
  *
  * Algorithm Overview:
  * 1. Range Reduction: Reduces input x to r where |r| < π/16
  *    - Computes n = round(x * 16/π) and r = x - n*π/16
- *    - Uses multi-precision arithmetic (3 parts of π/16) for accuracy
+ *    - Uses multi-word π/16 (3 words with FMA, 4 without) for accuracy
  *
  * 2. Table Lookup: Retrieves precomputed sin(n*π/16) and cos(n*π/16)
- *    - Includes high and low precision parts for cos values
+ *    - Stores high and low parts for both sin and cos (lows packed 32:32)
  *
  * 3. Polynomial Approximation: Computes sin(r) and cos(r)
  *    - sin(r) ≈ r * (1 + r²*P_sin(r²)) where P_sin is a minimax polynomial
@@ -142,17 +143,16 @@
   // Note: cos_lo and sin_lo are packed together (32 bits each) to save memory.
   // cos_lo can be used as-is since it's in the upper bits, sin_lo needs
   // extraction. The precision loss is negligible for the final result.
-  // see data/lut-inl.h.sol for the table generation code.
+  // see data/kpi16-inl.h.sol for the table generation code.
   V sin_lo = BitCast(d, ShiftLeft<32>(BitCast(du, cos_lo)));
 
   // Step 3: Multi-precision computation of remainder r
   // r = x - n*(π/16)_high
+  // Without native FMA the non-tail parts carry 27/25/29 bits, so the two
+  // leading products n*part are exact for |n| < ceil(2^24*16/π) = 85445660.
+  // part2 and the tail round, but the idioms below capture that rounding.
   constexpr auto kPiDiv16Prec29 = data::kPiDiv16Prec29<kNativeFMA>;
   V r_hi = NegMulAdd(n, Set(d, kPiDiv16Prec29[0]), x);
-  if constexpr (!kNativeFMA) {
-    // For F64, we need to handle the low precision part separately
-    r_hi = NegMulAdd(n, Set(d, kPiDiv16Prec29[3]), r_hi);
-  }
   const V pi16_med = Set(d, kPiDiv16Prec29[1]);
   const V pi16_lo = Set(d, kPiDiv16Prec29[2]);
   V r_med = NegMulAdd(n, pi16_med, r_hi);
@@ -162,6 +162,15 @@
   V term = NegMulAdd(pi16_med, n, Sub(r_hi, r_med));
   V r_lo = MulAdd(pi16_lo, n, Sub(r, r_med));
   r_lo = Sub(term, r_lo);
+  if constexpr (!kNativeFMA) {
+    // Fourth piece (48 bits at 2^-91). Reusing the one rounded product in both
+    // the subtraction and its error capture keeps (r, r_lo) an exact
+    // double-double even when r is tiny near k*π/2.
+    const V tail_prod = Mul(n, Set(d, kPiDiv16Prec29[3]));
+    const V r_prev = r;
+    r = Sub(r_prev, tail_prod);
+    r_lo = Add(r_lo, Sub(Sub(r_prev, r), tail_prod));
+  }
 
   // Step 4: Polynomial approximation
   V r2 = Mul(r, r);
@@ -210,13 +219,23 @@
   if constexpr (OP == Operation::kCos) {
     // Cosine reconstruction: cos_table - sin_table*remainder
     // Equivalent to: cos(a)*cos(r) - sin(a)*sin(r) but more efficient
-    V res_hi = NegMulAdd(r, sin_hi, cos_hi);  // cos_hi - r*sin_hi
-
-    // This captures the precision lost in the main computation
-    V r_sin_hi = Sub(cos_hi, res_hi);  // Extract high part of multiplication
-
-    // Handles rounding errors and adds sin_low contribution
-    V r_sin_low = MulSub(r, sin_hi, r_sin_hi);  // Compute multiplication error
+    V res_hi, r_sin_low;
+    if constexpr (kNativeFMA) {
+      res_hi = NegMulAdd(r, sin_hi, cos_hi);  // cos_hi - r*sin_hi
+
+      // This captures the precision lost in the main computation
+      V r_sin_hi = Sub(cos_hi, res_hi);  // Extract high part of multiplication
+
+      // Handles rounding errors and adds the low-part contribution
+      r_sin_low = MulSub(r, sin_hi, r_sin_hi);  // Compute multiplication error
+    } else {
+      // Dekker split stands in for the FMA idiom
+      V r_sin, r_sin_rest;
+      SplitMul(r, sin_hi, r_sin, r_sin_rest);
+      res_hi = Sub(cos_hi, r_sin);
+      V r_sin_hi = Sub(cos_hi, res_hi);
+      r_sin_low = Add(Sub(r_sin, r_sin_hi), r_sin_rest);
+    }
     V sin_low_corr = MulAdd(r, sin_lo, r_sin_low);  // Add sin_low term
 
     // This is used to apply the low-precision remainder correction
@@ -245,13 +264,23 @@
   } else {
     // Sine reconstruction: sin_table + cos_table*remainder
     // Equivalent to: sin(a)*cos(r) + cos(a)*sin(r) but more efficient
-    V res_hi = MulAdd(r, cos_hi, sin_hi);  // sin_hi + r*cos_hi
-
-    // This captures the precision lost in the main computation
-    V r_cos_hi = Sub(res_hi, sin_hi);  // Extract high part of multiplication
-
-    // Handles rounding errors and adds cos_low contribution
-    V r_cos_low = MulSub(r, cos_hi, r_cos_hi);  // Compute multiplication error
+    V res_hi, r_cos_low;
+    if constexpr (kNativeFMA) {
+      res_hi = MulAdd(r, cos_hi, sin_hi);  // sin_hi + r*cos_hi
+
+      // This captures the precision lost in the main computation
+      V r_cos_hi = Sub(res_hi, sin_hi);  // Extract high part of multiplication
+
+      // Handles rounding errors and adds the low-part contribution
+      r_cos_low = MulSub(r, cos_hi, r_cos_hi);  // Compute multiplication error
+    } else {
+      // Dekker split stands in for the FMA idiom
+      V r_cos, r_cos_rest;
+      SplitMul(r, cos_hi, r_cos, r_cos_rest);
+      res_hi = Add(sin_hi, r_cos);
+      V r_cos_hi = Sub(res_hi, sin_hi);
+      r_cos_low = Add(Sub(r_cos, r_cos_hi), r_cos_rest);
+    }
     V cos_low_corr = MulAdd(r, cos_lo, r_cos_low);  // Add cos_low term
 
     // Intermediate term for r_low correction: cos_table - sin_table*r

npsr/trig/inl.h

@@ -3,9 +3,10 @@
 // configurable precision, special case handling, and algorithm selection
 //
 // The implementation automatically selects between three algorithms:
-// 1. Low precision: ~1-4 ULP error, fastest
+// 1. Low precision: ~2-3 ULP error, fastest
 // 2. High precision: ~1 ULP error, moderate speed
-// 3. Extended precision: Exact for |x| > 2^24 (float) or 2^53 (double)
+// 3. Extended precision: Payne-Hanek reduction for huge |x|
+//    (> 10^4 for float, > 2^24 for double)
 
 #if defined(NPSR_TRIG_INL_H_) == defined(HWY_TARGET_TOGGLE)  // NOLINT
 #ifdef NPSR_TRIG_INL_H_
@@ -22,7 +23,7 @@
 
 HWY_BEFORE_NAMESPACE();
 
 namespace npsr::HWY_NAMESPACE::trig {
 /**
  * @brief Unified sine/cosine implementation with configurable precision
  *
@@ -46,7 +47,8 @@
  *
  * Thresholds for extended precision:
  * - Float: |x| > 10,000 (empirically chosen for accuracy)
- * - Double: |x| > 2^24 (16,777,216 - where 53-bit mantissa loses precision)
+ * - Double: |x| > 2^24 (16,777,216 - beyond this the reduction's n*π products
+ *   stop being exactly representable)
  */
 template <Operation OP, typename Prec, typename V>
 NPSR_INTRIN V Trig(Prec &prec, V x) {
@@ -56,8 +58,8 @@
   V ret;
   // Step 1: Select base algorithm based on accuracy requirements
   if constexpr (Prec::kLowAccuracy) {
-    // Low precision: Cody-Waite reduction with degree-9 polynomial
-    // Error: ~2 ULP and 3~ for non-fma
+    // Low precision: Cody-Waite reduction, degree-9 (f32) / degree-15 (f64)
+    // Error: ~2-3 ULP
     ret = Low<OP>(x);
   } else {
     // High precision: π/16 reduction with table lookup + polynomial
@@ -78,11 +80,16 @@
   // For |x| > threshold, standard algorithms lose precision due to
   // catastrophic cancellation in x - n*π reduction
   if constexpr (Prec::kLargeArgument) {
-    // Thresholds chosen based on when standard reduction loses accuracy:
-    // - Float: 10,000 is conservative but ensures < 1 ULP error
-    // - Double: 2^24 is where mantissa can't represent x and x-2π distinctly
+    // Thresholds mark where each path's n*π products stop being exact:
+    // - Float: 10,000, conservative but keeps the widened path < 1 ULP.
+    // - Double: 2^24, past which |round(x*16/π)| overflows the π/16 split.
+    // - Double low-accuracy cos: holds until trunc(|x|/π) > 8388606; the
+    //   constant is the largest double below that boundary (~8388607π).
+    constexpr bool kIsLowCos =
+        Prec::kLowAccuracy && OP == trig::Operation::kCos;
+    constexpr double kLargeF64 = kIsLowCos ? 0x1.921fb2200366fp+24 : 16777216.0;
     auto has_large_arg =
-        And(Gt(Abs(x), Set(d, kIsSingle ? 10000.0f : 16777216.0)), is_finite);
+        And(Gt(Abs(x), Set(d, kIsSingle ? 10000.0f : kLargeF64)), is_finite);
 
     // Extended precision is expensive, only use when necessary
     if (HWY_UNLIKELY(!AllFalse(d, has_large_arg))) {
@@ -101,7 +108,7 @@
 }  // namespace npsr::HWY_NAMESPACE::trig
 
 // Public API in the main npsr namespace
 namespace npsr::HWY_NAMESPACE {
 
 /**
  * @brief Compute sine of vector elements with configurable precision

npsr/trig/low-inl.h

@@ -10,7 +10,7 @@
 
 HWY_BEFORE_NAMESPACE();
 
 namespace npsr::HWY_NAMESPACE::trig {
 
 enum class Operation { kSin = 0, kCos = 1 };
 
@@ -56,13 +56,17 @@
   using namespace hn;
 
   const DFromV<V> d;
-  const V c15 = Set(d, -0x1.9f1517e9f65fp-41);
-  const V c13 = Set(d, 0x1.60e6bee01d83ep-33);
-  const V c11 = Set(d, -0x1.ae6355aaa4a53p-26);
-  const V c9 = Set(d, 0x1.71de3806add1ap-19);
-  const V c7 = Set(d, -0x1.a01a019a659ddp-13);
-  const V c5 = Set(d, 0x1.111111110a573p-7);
-  const V c3 = Set(d, -0x1.55555555554a8p-3);
+  // Both ops fit sin(r) on the same interval, so the two coefficient sets
+  // differ only in trailing bits; each matches its SVML counterpart
+  // (__svml_sin_d_la / __svml_cos_d_la) bit-for-bit.
+  constexpr bool kCos = OP == Operation::kCos;
+  const V c15 = Set(d, kCos ? -0x1.9f0d60811aac8p-41 : -0x1.9f1517e9f65fp-41);
+  const V c13 = Set(d, kCos ? 0x1.60e6857a2f220p-33 : 0x1.60e6bee01d83ep-33);
+  const V c11 = Set(d, kCos ? -0x1.ae63546002231p-26 : -0x1.ae6355aaa4a53p-26);
+  const V c9 = Set(d, kCos ? 0x1.71de38030fea0p-19 : 0x1.71de3806add1ap-19);
+  const V c7 = Set(d, kCos ? -0x1.a01a019a5b87bp-13 : -0x1.a01a019a659ddp-13);
+  const V c5 = Set(d, kCos ? 0x1.111111110a4a8p-7 : 0x1.111111110a573p-7);
+  const V c3 = Set(d, kCos ? -0x1.55555555554a7p-3 : -0x1.55555555554a8p-3);
   V poly = MulAdd(c15, r2, c13);
   poly = MulAdd(r2, poly, c11);
   poly = MulAdd(r2, poly, c9);
@@ -110,26 +114,28 @@
     //   N' = N - 0.5
     n = Sub(n, Set(d, static_cast<T>(0.5)));
   }
-  // Use Cody-Waite method with triple-precision PI
+  // Cody-Waite reduction with multi-word π (3 words with FMA, 4 without)
   constexpr auto kPi = data::kPi<T, kNativeFMA>;
-
   V r = NegMulAdd(n, Set(d, kPi[0]), x_abs);
   r = NegMulAdd(n, Set(d, kPi[1]), r);
   V r_lo = NegMulAdd(n, Set(d, kPi[2]), r);
+
   if constexpr (!kNativeFMA) {
-    if (!kIsSingle) {
-      r = r_lo;
-    }
     r_lo = NegMulAdd(n, Set(d, kPi[3]), r_lo);
   }
-
-  if constexpr (kIsSingle) {
+  if constexpr (kIsSingle || !kNativeFMA) {
+    // The polynomial needs the fully reduced value; r still owes the last
+    // π word(s).
     r = r_lo;
   }
+
   V r2 = Mul(r, r);
   V poly = PolyLow<OP>(r, r2);
 
   if constexpr (!kIsSingle) {
+    // Non-FMA has r == r_lo, giving the plain r + r³·poly. With FMA this is
+    // r_lo·(1 + r²·poly) ≈ sin(r_lo): r and r_lo differ by n·π[2] (~2^-83)
+    // and sin(r)/r varies slowly, so the error stays sub-ULP.
     V r2_corr = Mul(r2, r_lo);
     poly = MulAdd(r2_corr, poly, r_lo);
   }