ENH: improve f32 sin/cos accuracy

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,15 @@
 #define NPSR_LUT_INL_H_
 #endif
 
+#include <limits>
 #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, };
@@ -19,7 +19,7 @@
 template <bool FMA> inline constexpr double kPiPrec35[] = {
 0x1.921fb5444p1, 0x1.68c234c4c6628p-38, };
 template <> inline constexpr double kPiPrec35<false>[] = {
-0x1.921fb6p1, -0x1.777a5cp-24, -0x1.ee59dap-49, };
+0x1.921fb544p1, 0x1.0b4611a8p-33, -0x1.d9cceba3f91f2p-65, };
 
 template <typename T> inline constexpr char kPiMul2[] = {};
 template <> inline constexpr float kPiMul2<float>[] = {
@@ -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/data/constants.h.sol

@@ -44,12 +44,20 @@ Append(
   // rounds ~2^-105, harmless even for tiny reduced arguments near k*π/2
   "",
 
-  // Special 35-bit precision π for specific algorithms
+  // Special 35-bit precision π for the f32 High path, which reduces in double.
+  // FMA: the fused subtractions never need n*πᵢ to be representable, so two
+  // words carry the whole budget and π is held to 35+53 = 88 bits.
+  //
+  // non-FMA: NegMulAdd decays to `x - fl(n*πᵢ)`, so each n*πᵢ has to be exact.
+  // The quotient is bounded by the f32 magic-round in High(), |n| < 2^22, hence
+  // 53-22 = 31-bit heads. Two of those only pin π to 62 bits, so a 53-bit tail
+  // follows; its product rounds at ~2^-92, just under the FMA split's 2^-90
+  // residual, which is what keeps both variants on the same accuracy curve.
   "template <bool FMA> inline constexpr double kPiPrec35[] = " @
   KArray_(Float64, pi, [|RN, 35|], [|RD, 53|]),
-  
+
   "template <> inline constexpr double kPiPrec35<false>[] = " @
-  KArray_(Float64, pi, [|RN, 24, 24, 24|]),
+  KArray_(Float64, pi, [|RN, 31, 31|], [|RN, 53|]),
   "",
 
   // 2π constants for angle wrapping

npsr/trig/data/data.h

@@ -1,3 +1,3 @@
 // Aggregates the trig data headers. Hand-written: it only forwards includes,
 // so there is nothing for Sollya to generate.
 //
@@ -6,8 +6,9 @@
 // hwy/foreach_target.h. An include-once guard would suppress all but the first
 // target pass. The include-once children (constants/approx/reduction) carry
 // their own guards and no-op on re-entry.
 #include "npsr/lut-inl.h"
 #include "npsr/trig/data/constants.h"
+#include "npsr/trig/data/polyf32.h"
 #include "npsr/trig/data/kpi16-inl.h"
 #include "npsr/trig/data/approx.h"
 #include "npsr/trig/data/reduction.h"

npsr/trig/data/polyf32.h

@@ -0,0 +1,32 @@
+// Auto-generated by npsr/trig/data/polyf32.h.sol
+// Use `spin sollya -f` to force regeneration
+#ifndef NPSR_TRIG_DATA_POLYF32_H
+#define NPSR_TRIG_DATA_POLYF32_H
+
+namespace npsr::trig::data {
+// Degree-9 odd minimax sin(x) = x + x^3*P(x^2) on [2^-20, pi/2].
+// Coefficients {c3, c5, c7, c9} for the f32 PolyLow in npsr/trig/low-inl.h.
+// See poly.h.sol for the recipe, pins and measured worst cases.
+template <bool FMA> inline constexpr float kSinPolyLowF32[] = {
+-0x1.555544p-3f, 0x1.110d94p-7f, -0x1.9f3644p-13f, 0x1.5a9eap-19f, 
+};
+template <bool FMA> inline constexpr float kCosPolyLowF32[] = {
+-0x1.555548p-3f, 0x1.110e34p-7f, -0x1.9f527cp-13f, 0x1.5c1b4p-19f, 
+};
+template <> inline constexpr float kSinPolyLowF32<false>[] = {
+-0x1.55554ap-3f, 0x1.110e88p-7f, -0x1.9f621ap-13f, 0x1.5cfce2p-19f, 
+};
+template <> inline constexpr float kCosPolyLowF32<false>[] = {
+-0x1.55554ap-3f, 0x1.110e9p-7f, -0x1.9f653ep-13f, 0x1.5d3afep-19f, 
+};
+
+// Degree-9 odd minimax sin(r) on [2^-30, pi/2 + 0.012] for the f32 High path
+// in trig/high-inl.h; c1 is fitted, not pinned to 1. See poly.h.sol.
+// Relative error 9.64818894863128662109375e-2 ulp(f32). Coefficients {c1, c3, c5, c7, c9}.
+inline constexpr double kSinPolyHighF32[] = {
+0x1.ffffffce99ecep-1, -0x1.5555471207a17p-3, 0x1.110e5bbaaff44p-7, 
+-0x1.9f5e996744e9ep-13, 0x1.5ce7e229f1f1p-19, };
+
+} // namespace npsr::trig::data
+
+#endif // NPSR_TRIG_DATA_POLYF32_H

npsr/trig/data/polyf32.h.sol

@@ -0,0 +1,113 @@
+// f32 Low trig coefficients: sin(x) ~= x + x^3*P(x^2), degree 9, odd.
+// 333/346/368 are recoverable Remez warnings the constrained refits below
+// always trip; the driver would treat them as errors.
+suppressmessage(346, 333, 368);
+
+I = [0x1p-20; pi/2];
+vFree = fpminimax(sin(x), [|3, 5, 7, 9|], [|single...|], I, relative, floating, x);
+
+// Relative error in ULP at 1 (= 2^-24); guards regeneration against a worse fit.
+procedure ErrULP(pP) { return dirtyinfnorm(pP(x)/sin(x) - 1, I) * 0x1p24; };
+Assert(ErrULP(vFree) < 0.12, "f32 Low free fit error regressed: " @ ErrULP(vFree));
+
+// One f32 ULP at a coefficient's scale. No coefficient sits within a step of
+// a power of two, so the constant step is exact.
+procedure Ulp(pC) { return 2^(floor(log2(abs(pC))) - 23); };
+
+// Pin c3 (then c5, then optionally c7) an integer number of ULPs off the free
+// fit and refit the rest under the pin. Pin indices come from exhaustive
+// measurement of the compiled kernel over [0, 1e4); no analytic criterion
+// sees the eval DAG's per-operation rounding.
+procedure Stage2(pK3, pJ5) {
+  var vC3, vP5, vC5, vP7;
+  vC3 = round(coeff(vFree, 3) + pK3 * Ulp(coeff(vFree, 3)), single, RN);
+  vP5 = fpminimax(sin(x), [|5, 7, 9|], [|single...|], I, relative, floating,
+                  x + vC3 * x^3);
+  vC5 = round(coeff(vP5, 5) + pJ5 * Ulp(coeff(vP5, 5)), single, RN);
+  vP7 = fpminimax(sin(x), [|7, 9|], [|single...|], I, relative, floating,
+                  x + vC3 * x^3 + vC5 * x^5);
+  return [| vC3, vC5, vP7 |];
+};
+
+procedure CheckPinned(pSet, pTag) {
+  Assert(ErrULP(pSet[0]*x^3 + pSet[1]*x^5 + pSet[2]*x^7 + pSet[3]*x^9 + x) < 1.0,
+         "pinned set " @ pTag @ " analytic error drifted");
+  return pSet;
+};
+
+// c7/c9 come jointly from the stage-two fit.
+procedure Pinned(pK3, pJ5) {
+  var vS;
+  vS = Stage2(pK3, pJ5);
+  return CheckPinned([| vS[0], vS[1], coeff(vS[2], 7), coeff(vS[2], 9) |],
+                     "(" @ pK3 @ "," @ pJ5 @ ")");
+};
+
+// Additionally pin c7 and refit c9 alone. Only FMA cos gains from it.
+procedure Pinned3(pK3, pJ5, pL7) {
+  var vS, vC7, vP9;
+  vS = Stage2(pK3, pJ5);
+  vC7 = round(coeff(vS[2], 7) + pL7 * Ulp(coeff(vS[2], 7)), single, RN);
+  vP9 = fpminimax(sin(x), [|9|], [|single...|], I, relative, floating,
+                  x + vS[0] * x^3 + vS[1] * x^5 + vC7 * x^7);
+  return CheckPinned([| vS[0], vS[1], vC7, coeff(vP9, 9) |],
+                     "(" @ pK3 @ "," @ pJ5 @ "," @ pL7 @ ")");
+};
+
+// Exhaustive [0, 1e4) worst case, pinned vs free fit, in ULP:
+//   FMA:     sin 1.886 (free 2.191)   cos 1.815 (free 2.110)
+//   non-FMA: sin 1.912 (free 1.939)   cos 1.948 (free 1.966)
+// Non-FMA is eval-rounding bound (~1.9 of the ~1.95 total), so pins only
+// shift which argument realizes the worst alignment; the whole (k3, j5, l7)
+// grid stays within 0.02 of the sets below.
+vFmaSin = Pinned(4, 1);
+vFmaCos = Pinned3(2, 1, 8);
+vNoFmaSin = Pinned(1, 1);
+vNoFmaCos = Pinned(1, 5);
+
+Append(
+  "// Degree-9 odd minimax sin(x) = x + x^3*P(x^2) on [2^-20, pi/2].",
+  "// Coefficients {c3, c5, c7, c9} for the f32 PolyLow in npsr/trig/low-inl.h.",
+  "// See poly.h.sol for the recipe, pins and measured worst cases.",
+  "template <bool FMA> inline constexpr float kSinPolyLowF32[] = " @ CArrayT(Float32, vFmaSin, 4) @ ";",
+  "template <bool FMA> inline constexpr float kCosPolyLowF32[] = " @ CArrayT(Float32, vFmaCos, 4) @ ";",
+  "template <> inline constexpr float kSinPolyLowF32<false>[] = " @ CArrayT(Float32, vNoFmaSin, 4) @ ";",
+  "template <> inline constexpr float kCosPolyLowF32<false>[] = " @ CArrayT(Float32, vNoFmaCos, 4) @ ";",
+  ""
+);
+
+// f32 High coefficients, held and evaluated in double by npsr/trig/high-inl.h.
+//
+// c1 is fitted, not pinned to 1: MulAdd(r2, poly, c1) + Mul(poly, r) closes in
+// the same two ops as the pinned x + x^3*P(x^2), so the free parameter is free
+// -- 0.1099 -> 0.0965 ulp. No pinning stage either: eval is in double (~2^-52),
+// so unlike the single Low sets above there is no rounding for pins to steer.
+//
+// The interval tracks the cutover in npsr/trig/inl.h and must move with it: n
+// comes from a f32 magic-round, so |r| overshoots pi/2 by up to |x|*1.7e-7.
+// 0.012 is the exhaustive ceiling below both cutovers, and widening costs
+// accuracy fast: 0.0892 ulp at pi/2, 0.0965 here, 0.1073 at 0.028.
+IHigh = [0x1p-30; pi/2 + 0.012];
+vHigh = fpminimax(sin(x), [|1, 3, 5, 7, 9|], [|double...|], IHigh, relative, floating);
+
+// Leaves High<> at 0.5 (final f32 rounding) + 0.097, the f32 ceiling: the
+// Extended kernel above the cutover measures 0.5031, so this fit -- not the
+// Payne-Hanek reduction -- is what keeps sin/cos off correct rounding.
+errHigh = dirtyinfnorm(vHigh(x)/sin(x) - 1, IHigh) * 0x1p24;
+Assert(errHigh < 0.0970, "f32 High fit error regressed: " @ errHigh);
+// The driver leaves display in hexadecimal, which is unreadable for an error.
+display = decimal!;
+errHighStr = "" @ round(errHigh, 24, RN);
+display = hexadecimal!;
+
+Append(
+  "// Degree-9 odd minimax sin(r) on [2^-30, pi/2 + 0.012] for the f32 High path",
+  "// in trig/high-inl.h; c1 is fitted, not pinned to 1. See poly.h.sol.",
+  "// Relative error " @ errHighStr @ " ulp(f32). Coefficients {c1, c3, c5, c7, c9}.",
+  "inline constexpr double kSinPolyHighF32[] = " @
+  CArrayT(Float64, [| coeff(vHigh, 1), coeff(vHigh, 3), coeff(vHigh, 5),
+                      coeff(vHigh, 7), coeff(vHigh, 9) |], 3) @ ";",
+  ""
+);
+
+WriteCPPHeader("npsr::trig::data");

npsr/trig/extended-inl.h

@@ -10,7 +10,7 @@
 #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)