BUG: Fix f64 trig accuracy issues

- Match SVML rounding DAG in large-argument reconstruction
- Give low-accuracy f64 cosine its own minimax fit (Match SVML)
- Restore sub-ULP accuracy on non-FMA targets
- Mark dead x-sign term in High sign DAG for removal

Author: seiko2plus

npsr/lut-inl.h

@@ -6,6 +6,7 @@
 #endif
 
 #include <tuple>
+#include <limits>
 
 #include "npsr/hwy.h"
 

npsr/trig/data/approx.h

@@ -30,7 +30,7 @@ template <> inline constexpr double kPiMul2<double>[] = {
 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/data/constants.h

@@ -233,6 +233,16 @@ NPSR_INTRIN V Extended(V x) {
     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,31 @@ NPSR_INTRIN V Extended(V x) {
   // =============================================================================
   // 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: near k*pi/2, res_hi is
+    // ~deriv_hi*r, so the product rounding lost by a decayed MulAdd costs
+    // up to 0.5 ulp of the result. (sigma*r in the chains below is fine
+    // decayed: sigma ~ 2^-53*deriv makes its product error second-order.)
+    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 +320,22 @@ NPSR_INTRIN V Extended(V x) {
   }
   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/extended-inl.h

@@ -147,12 +147,13 @@ NPSR_INTRIN V High(V x) {
 
   // 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/27/29 bits so every
+  // n*part product is exact for |n| <= 85445660 = round(2^24 * 16/π); each
+  // subtraction's rounding error is then captured exactly below, and 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,6 +163,15 @@ NPSR_INTRIN V High(V x) {
   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 (53 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);
@@ -210,13 +220,23 @@ NPSR_INTRIN V High(V x) {
   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 +265,23 @@ NPSR_INTRIN V High(V x) {
   } 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,6 +312,9 @@ NPSR_INTRIN V High(V x) {
   // 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
   VU combined = Xor(BitCast(du, n_biased), ShiftLeft<4>(x_sign_int));

npsr/trig/high-inl.h

@@ -49,7 +49,7 @@ namespace npsr::HWY_NAMESPACE::trig {
  * - Double: |x| > 2^24 (16,777,216 - where 53-bit mantissa loses precision)
  */
 template <Operation OP, typename Prec, typename V>
-NPSR_INTRIN V Trig(Prec &prec, V x) {
+NPSR_INTRIN V Trig(Prec& prec, V x) {
   using namespace hwy::HWY_NAMESPACE;
   constexpr bool kIsSingle = std::is_same_v<TFromV<V>, float>;
   const DFromV<V> d;
@@ -81,8 +81,14 @@ NPSR_INTRIN V Trig(Prec &prec, V x) {
     // 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
+    // - Double low-accuracy cosine:  stays on its fast path
+    //   until trunc(|x|/π) > 8388606; this 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))) {
@@ -121,7 +127,7 @@ namespace npsr::HWY_NAMESPACE {
  * ```
  */
 template <typename Prec, typename V>
-NPSR_INTRIN V Sin(Prec &prec, V x) {
+NPSR_INTRIN V Sin(Prec& prec, V x) {
   return trig::Trig<trig::Operation::kSin>(prec, x);
 }
 
@@ -143,7 +149,7 @@ NPSR_INTRIN V Sin(Prec &prec, V x) {
  * ```
  */
 template <typename Prec, typename V>
-NPSR_INTRIN V Cos(Prec &prec, V x) {
+NPSR_INTRIN V Cos(Prec& prec, V x) {
   return trig::Trig<trig::Operation::kCos>(prec, x);
 }
 

npsr/trig/inl.h

@@ -14,6 +14,30 @@ namespace npsr::HWY_NAMESPACE::trig {
 
 enum class Operation { kSin = 0, kCos = 1 };
 
+/**
+ * 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));
+}
+
 template <Operation OP, typename V, HWY_IF_F32(TFromV<V>)>
 NPSR_INTRIN V PolyLow(V r, V r2) {
   using namespace hn;
@@ -54,15 +78,17 @@ NPSR_INTRIN V PolyLow(V r, V r2) {
 template <Operation OP, typename V, HWY_IF_F64(TFromV<V>)>
 NPSR_INTRIN V PolyLow(V r, V r2) {
   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);
+  // Sine and cosine use distinct minimax fits of sin(r): the cosine path
+  // reduces around (N + 0.5)*pi, so its coefficient set differs 
+  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);
@@ -117,10 +143,14 @@ NPSR_INTRIN V Low(V x) {
   r = NegMulAdd(n, Set(d, kPi[1]), r);
   V r_lo = NegMulAdd(n, Set(d, kPi[2]), r);
   if constexpr (!kNativeFMA) {
-    if (!kIsSingle) {
+    r_lo = NegMulAdd(n, Set(d, kPi[3]), r_lo);
+    if constexpr (!kIsSingle) {
+      // The polynomial must consume the fully reduced value: the partial
+      // chain still carries ~n*2^-49 of unapplied π tail, far too coarse
+      // for r² and its correction terms. (The f32 path already reads the
+      // chain's end below.)
       r = r_lo;
     }
-    r_lo = NegMulAdd(n, Set(d, kPi[3]), r_lo);
   }
 
   if constexpr (kIsSingle) {