ENH: improve f32 sin/cos accuracy

seiko2plus · seiko2plus · commit 25c78204639f · 2026-07-28T03:27:18.000+03:00
The new improvements provide better accuracy for both non-FMA and FMA
targets, so we can achieve &lt; 2ulp for the low accuracy path and &lt; 0.6 for
the high accuracy path without affecting performance, except for a few ops on
non-FMA targets, which was necessary to match FMA target precision.
The improvements also include increasing the fast path ranges for FMA targets,
and it has been tested against an exhaustive test covering all float32 ranges.

This was achieved by refitting the f32 sin/cos polynomials as generated minimax sets (new
data/polyf32.h.pol) instead of SVML hardcoded coefficients, split per op and per
FMA mode. For non-FMA:
- High&lt;&gt; gets a third pi word with 31-bit heads so n*pi_i stays
- exact; FMA keeps the two-word 35+53-bit split.
- Low&lt;&gt; recovers the pi[2] rounding error via Fast2Sum on non-FMA,
  the same approach applies to float64 and buys ~0.3 ulp.

Also, for all targets:
- Low&lt;&gt; cos(x) == 1 is answered for all |x| below a per-type bound, not just
x_trans == pi/2, which fixes tiny tests.
- kLargeF32 replaces the flat 10000 with per-path measured thresholds.
- Sollya prec 512 -&gt; 1024, which was needed for generating the float32-widening double poly.
diff --git a/npsr/lut-inl.h b/npsr/lut-inl.h
@@ -7,6 +7,7 @@
 
 #include <limits>
 #include <tuple>
+#include <limits>
 
 #include "npsr/hwy.h"
 
diff --git a/npsr/trig/data/constants.h b/npsr/trig/data/constants.h
@@ -19,7 +19,7 @@ template <> inline constexpr double kPi<double, false>[] = {
 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>[] = {
diff --git a/npsr/trig/data/constants.h.sol b/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
diff --git a/npsr/trig/data/data.h b/npsr/trig/data/data.h
@@ -8,6 +8,7 @@
 // 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"
diff --git a/npsr/trig/data/polyf32.h b/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
diff --git a/npsr/trig/data/polyf32.h.sol b/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");
diff --git a/npsr/trig/high-inl.h b/npsr/trig/high-inl.h
@@ -59,21 +59,27 @@ NPSR_INTRIN V High(V x) {
   }
   auto WideCal = [](const VW &nh, const VW &xh_abs) -> VW {
     const DFromV<VW> dw;
-    constexpr auto kPiPrec35 = data::kPiPrec35<true>;
+    // Cody-Waite r = x - n*π, one word per NegMulAdd. Without FMA the product
+    // rounds, so that split spends a third word on 31-bit heads to keep n*πᵢ
+    // exact; see constants.h.sol.
+    constexpr auto kPiPrec35 = data::kPiPrec35<kNativeFMA>;
     VW r = NegMulAdd(nh, Set(dw, kPiPrec35[0]), xh_abs);
     r = NegMulAdd(nh, Set(dw, kPiPrec35[1]), r);
+    if constexpr (!kNativeFMA) {
+      r = NegMulAdd(nh, Set(dw, kPiPrec35[2]), r);
+    }
     VW r2 = Mul(r, r);
 
-    // Polynomial coefficients for sin(r) approximation on [-π/2, π/2]
-    const VW c9 = Set(dw, 0x1.5dbdf0e4c7deep-19);
-    const VW c7 = Set(dw, -0x1.9f6ffeea73463p-13);
-    const VW c5 = Set(dw, 0x1.110ed3804ca96p-7);
-    const VW c3 = Set(dw, -0x1.55554bc836587p-3);
-    VW poly = MulAdd(c9, r2, c7);
-    poly = MulAdd(r2, poly, c5);
-    poly = MulAdd(r2, poly, c3);
-    poly = Mul(poly, r2);
-    poly = MulAdd(r, poly, r);
+    // Degree-9 odd minimax for sin(r); the fit is this path's whole error
+    // budget (0.5 + 0.0965 ULP). c1 is fitted, not pinned to 1 -- same two
+    // closing ops as r + r^3*P(r^2), one more free parameter. Within 2^-27.4
+    // of 1, so sin(r) still demotes to r bit-exactly. See data/polyf32.h.sol.
+    constexpr auto kPoly = data::kSinPolyHighF32;
+    VW poly = MulAdd(Set(dw, kPoly[4]), r2, Set(dw, kPoly[3]));
+    poly = MulAdd(r2, poly, Set(dw, kPoly[2]));
+    poly = MulAdd(r2, poly, Set(dw, kPoly[1]));
+    poly = MulAdd(r2, poly, Set(dw, kPoly[0]));
+    poly = Mul(poly, r);
     return poly;
   };
 
diff --git a/npsr/trig/inl.h b/npsr/trig/inl.h
@@ -1,13 +1,12 @@
-// Main trigonometric function dispatcher for Highway SIMD library
 // This file provides the public API for sine and cosine functions with
-// configurable precision, special case handling, and algorithm selection
+// configurable precision, special case handling, and algorithm selection.
+//
+// The implementation automatically selects between:
+// 1. Low precision: < 2 ULP for f32 and < 3.5 ULP for f64
+// 2. High precision: < 0.6 ULP for f32, < 1 ULP for f64; moderate speed
+// 3. Extended precision: Payne-Hanek reduction for huge |x|, similar
+// accuracy to High precision but a slow path.
 //
-// The implementation automatically selects between three algorithms:
-// 1. Low precision: ~2-3 ULP error, fastest
-// 2. High precision: ~1 ULP error, moderate speed
-// 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_
 #undef NPSR_TRIG_INL_H_
@@ -42,11 +41,13 @@ namespace npsr::HWY_NAMESPACE::trig {
  *
  * Algorithm selection:
  * 1. If kLowAccuracy: Use Low<> (Cody-Waite with minimal polynomial)
- * 2. Otherwise: Use High<> (π/16 reduction with table lookup)
+ * 2. Otherwise: Use High<> (f32: π reduction in double; f64: π/16 reduction
+ *    with table lookup)
  * 3. If kLargeArgument and |x| > threshold: Override with Extended<>
  *
  * Thresholds for extended precision:
- * - Float: |x| > 10,000 (empirically chosen for accuracy)
+ * - Float: |x| > the largest argument the chosen path still meets its bound at
+ *   (differs per accuracy/FMA/op; see kLargeF32)
  * - Double: |x| > 2^24 (16,777,216 - beyond this the reduction's n*π products
  *   stop being exactly representable)
  */
@@ -58,12 +59,8 @@ NPSR_INTRIN V Trig(Prec &prec, V x) {
   V ret;
   // Step 1: Select base algorithm based on accuracy requirements
   if constexpr (Prec::kLowAccuracy) {
-    // 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
-    // Error: ~1 ULP
     ret = High<OP>(x);
   }
   // Step 2: Handle special cases (NaN, Inf) if enabled
@@ -80,16 +77,26 @@ NPSR_INTRIN V Trig(Prec &prec, V x) {
   // For |x| > threshold, standard algorithms lose precision due to
   // catastrophic cancellation in x - n*π reduction
   if constexpr (Prec::kLargeArgument) {
-    // Thresholds mark where each path's n*π products stop being exact:
-    // - Float: 10,000, conservative but keeps the widened path < 1 ULP.
+    // Thresholds mark where each path stops meeting its bound, i.e. where its
+    // range reduction breaks down, not where mantissa resolution runs out:
+    // - Float: each path keeps its own accuracy bound up to a different |x|;
+    //   every constant below is the round value under the exhaustively
+    //   measured, MPFR-confirmed first argument that exceeds it. The low paths
+    //   break far earlier without FMA, because there the extra π word is spent
+    //   on keeping n*πᵢ representable rather than on headroom.
     // - 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 bool kIsCos = OP == trig::Operation::kCos;
+    constexpr bool kIsLowCos = Prec::kLowAccuracy && kIsCos;
+    constexpr float kLargeF32 =
+        Prec::kLowAccuracy
+            ? (kNativeFMA ? (kIsCos ? 1600000.0f : 2200000.0f)
+                          : (kIsCos ? 12000.0f : 25000.0f))
+            : (kNativeFMA ? 200000.0f : 130000.0f);
     constexpr double kLargeF64 = kIsLowCos ? 0x1.921fb2200366fp+24 : 16777216.0;
     auto has_large_arg =
-        And(Gt(Abs(x), Set(d, kIsSingle ? 10000.0f : kLargeF64)), is_finite);
+        And(Gt(Abs(x), Set(d, kIsSingle ? kLargeF32 : kLargeF64)), is_finite);
 
     // Extended precision is expensive, only use when necessary
     if (HWY_UNLIKELY(!AllFalse(d, has_large_arg))) {
diff --git a/npsr/trig/low-inl.h b/npsr/trig/low-inl.h
diff --git a/tools/sollya/core.sol b/tools/sollya/core.sol
