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,30 @@
 using hn::VFromD;
 constexpr bool kNativeFMA = HWY_NATIVE_FMA != 0;
 
-inline HWY_ATTR void DummyToSuppressUnusedWarning() {}
+/**
+ * Substitute for the FMA product-error idiom on targets without native FMA.
+ *
+ * With FMA, `e = MulSub(a, b, Round(a*b))` recovers the rounding error of the
+ * product exactly. Without it, that expression only sees the already-rounded
+ * product and returns nothing useful. Instead, split both operands with a bit
+ * mask (26 significant bits in the head) so `head = a_hi*b_hi` is exact by
+ * construction, and return the dropped cross terms in `rest`:
+ *   a*b == head + rest, up to a negligible ~2^-105-scale rounding of `rest`.
+ */
+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 <tuple>
+#include <limits>
 
 #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

@@ -1,3 +1,3 @@
 #if defined(NPSR_TRIG_EXTENDED_INL_H_) == defined(HWY_TARGET_TOGGLE)  // NOLINT
 #ifdef NPSR_TRIG_EXTENDED_INL_H_
 #undef 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) {
@@ -41,7 +41,7 @@
   // providing ~96-bit (F32) or ~192-bit (F64) precision for (4/π) × 2^exp
   VU u_exponent = GetBiasedExponent(x);
   VI i_table_idx =
-      BitCast(di, Add(ShiftLeft<1>(u_exponent), u_exponent));  // × 2 + 1 = × 3
+      BitCast(di, Add(ShiftLeft<1>(u_exponent), u_exponent));  // 2·e + e = 3·e (3 words/entry)
 
   // Gather three parts of (4/π) × 2^exp from precomputed table
   // Generated by npsr/trig/data/reduction.h.sol
@@ -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,16 @@
     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) {
+    // F64 without FMA: `MulSub(pi2_hi, n, r)` sees only the already-rounded
+    // product and yields exactly zero, dropping the 0.5-ulp(r) rounding
+    // error of r. Near k*pi/2 the result is ~deriv*r, so that loss is worth
+    // 0.5 ulp of the final result. Recover it with the exact head product:
+    // head + rest == pi2_hi*n, and 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 +282,33 @@
   // =============================================================================
   // 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 FMA-idiom substitution as in Phase 7. res_hi = func_hi + sigma*r +
+    // deriv_hi*r must stay the faithfully-rounded head of the (res_hi, res_lo)
+    // pair: sigma is a signed power of two (see approx.h), so sigma*r is exact
+    // under a plain Mul and needs no capture, but deriv_hi*r would round under
+    // a decayed MulAdd and corrupt the head. SplitMul recovers that product
+    // rounding; the addition rounding is captured just below by Sub.
+    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 +322,22 @@
   }
   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);
+  // SVML pairs the correction terms differently per precision: F32 scales the
+  // sin term by corr (first-order cos at the reduced point) and accumulates
+  // the cos term onto the r_lo chain; F64 scales the sin term by deriv and
+  // accumulates the cos term onto it, adding the r_lo chain last.
+  V res_lo;
+  if constexpr (kIsSingle) {
+    V res_lo_sin = MulAdd(corr, sin_poly, res_lo0);
+    V res_lo_cos = MulAdd(func_hi, cos_poly, res_lo_rlo);
+    res_lo = Add(res_lo_cos, res_lo_sin);
+  } else {
+    V res_lo_sin = MulAdd(deriv, sin_poly, res_lo0);
+    V res_lo_cos = MulAdd(func_hi, cos_poly, res_lo_sin);
+    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) {
@@ -81,7 +81,9 @@
   VW poly_up = WideCal(PromoteUpperTo(dw, n), PromoteUpperTo(dw, x_abs));
 
   V poly = Combine(d, DemoteTo(dh, poly_up), DemoteTo(dh, poly_lo));
-  // Extract octant sign information from quotient and flip the sign bit
+  // sin(r) was reduced around n·π, and sin/cos flip sign every π, so the sign
+  // is the parity of n, carried in bit 0 of n_biased. Shift it into the sign
+  // position and XOR to flip the result on odd n.
   poly = Xor(poly,
              BitCast(d, ShiftLeft<sizeof(T) * 8 - 1>(BitCast(du, n_biased))));
   if constexpr (OP == Operation::kCos) {
@@ -95,15 +97,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 +145,20 @@
   // 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
+  // Both splits apply their parts in descending order through the same DAG.
+  // Without native FMA the non-tail parts carry 27/25/29 bits. The two
+  // leading products n*part are exact for |n| < 85445660 = ceil(2^24*16/π)
+  // (split and bound co-tuned so n*part0 just fits in 53 bits); part2 and
+  // the tail are rounded products, but their product rounding — like each
+  // subtraction's — is captured exactly by the idioms below, so the decayed
+  // mul+sub sequence loses nothing against the fused path.
   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,21 +168,31 @@
   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). Reuse the one rounded product in both
+    // the subtraction and its error capture so (r, r_lo) stays an exact
+    // double-double of the reduction 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);
 
-  // Minimax polynomial for (sin(r)/r - 1)
-  // sin(r)/r = 1 - r²/3! + r⁴/5! - r⁶/7! + ...
-  // This polynomial computes the terms after 1
+  // Minimax polynomial for (sin(r)/r - 1)/r²  ==  (sin(r) - r)/r³
+  // The Taylor form sin(r)/r - 1 = -r²/3! + r⁴/5! - ... fixes the leading
+  // coefficient at -1/3! = -1/6 (higher terms are minimax-adjusted); the
+  // stripped r² factor is reapplied downstream via MulAdd(r2, ...).
   V sin_poly = Set(d, 0x1.71c97d22a73ddp-19);
   sin_poly = MulAdd(sin_poly, r2, Set(d, -0x1.a01a00ed01edep-13));
   sin_poly = MulAdd(sin_poly, r2, Set(d, 0x1.111111110e99dp-7));
   sin_poly = MulAdd(sin_poly, r2, Set(d, -0x1.5555555555555p-3));
 
   // Minimax polynomial for (cos(r) - 1)/r²
-  // cos(r) = 1 - r²/2! + r⁴/4! - r⁶/6! + ...
-  // This polynomial computes (cos(r) - 1)/r²
+  // From cos(r) = 1 - r²/2! + r⁴/4! - ..., the leading term is ≈ -1/2! = -1/2
+  // (minimax-adjusted, so not exactly -1/2).
   V cos_poly = Set(d, 0x1.9ffd7d9d749bcp-16);
   cos_poly = MulAdd(cos_poly, r2, Set(d, -0x1.6c16c075d73f8p-10));
   cos_poly = MulAdd(cos_poly, r2, Set(d, 0x1.555555554e8d6p-5));
@@ -192,7 +208,7 @@
   //
   // Where for small r (|r| < π/16):
   // cos(r) ≈ 1 + r²*cos_poly
-  // sin(r) ≈ r*(1 + sin_poly) ≈ r + r*sin_poly
+  // sin(r) ≈ r + r³*sin_poly
   //
   // SVML's efficient linear approximation:
   // sin(a+r) ≈ sin(a) + cos(a)*r + polynomial_corrections
@@ -201,7 +217,7 @@
   // This is mathematically equivalent but computationally more efficient:
   // - Uses first-order linear terms directly: Sh + Ch*R, Ch - R*Sh
   // - Applies higher-order polynomial corrections separately
-  // - Fewer multiplications and better numerical stability
+  // - Better accuracy near the result's zeros (dominant term stays leading)
   //
   // Implementation follows SVML structure:
   // sin(n*π/16 + r) = sin_table + cos_table*remainder (+ corrections)
@@ -210,13 +226,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 via the FMA product-error idiom
+      r_sin_low = MulSub(r, sin_hi, r_sin_hi);  // Compute multiplication error
+    } else {
+      // Exact-by-construction head product replaces 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 +271,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 via the FMA product-error idiom
+      r_cos_low = MulSub(r, cos_hi, r_cos_hi);  // Compute multiplication error
+    } else {
+      // Exact-by-construction head product replaces 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
@@ -282,8 +318,11 @@
   // This unified approach works because:
   // - sin(x + π) = -sin(x)
   // - cos(x + π) = -cos(x)
+  // TODO(seiko2plus): x_sign_int and the Xor are dead code (always zero) and
+  // should be removed. The sign is bit 4 of n_biased alone; sin(-0.0) is
+  // handled in trig/inl.h under kSpecialCases.
   VU x_sign_int = ShiftLeft<63>(BitCast(du, x));
-  // XOR with quadrant info in n_biased
+  // Fold in the π-period parity (bit 4 of n_biased); x_sign_int is 0 (dead).
   VU combined = Xor(BitCast(du, n_biased), ShiftLeft<4>(x_sign_int));
   // Extract final sign
   VU sign = ShiftRight<4>(combined);