BUG: Regenerate trig LUT/constants from updated Sollya sources

Match SVML rounding DAG in the reduction LUT and use exact-product
constants for the non-FMA f64 path.

Author: seiko2plus

npsr/trig/data/approx.h.sol

@@ -1,79 +1,104 @@
-// Suppress rounding mode information messages and NaN warnings for boundary angles
+// Suppress rounding mode information messages for irrational table values
 suppressmessage(184, 185, 186); // suppress info no rounding, round-up, round-down
-suppressmessage(419); // expected nan when angle is at specific multiples causing derivative singularities
+suppressmessage(160);           // inequality decided by faithful evaluation
 
 // Generates a 4-element lookup table for fast trigonometric function approximation.
-// 
-// This procedure creates optimized lookup tables that store:
-// - Function values split into high/low precision components
-// - Derivative information with power-of-2 scaling for efficient interpolation
 //
 // The table format per angle is: [deriv, sigma, high, low]
 // where:
-//   deriv = actual_derivative - sigma (reduces storage requirements)
-//   sigma = 2^k, where k = ceil(log2(|actual_derivative|))
-//   high  = main part of function value
-//   low   = residual for extended precision
+//   sigma + deriv = f'(angle)  sigma = sign * 2^k nearest the derivative;
+//                              being a power of two, sigma*z is exact for any
+//                              small remainder z, confining the interpolation
+//                              rounding error to the tiny deriv*z term
+//   high  + low   = f(angle)   double-float split of the function value
+//
+// The values reproduce Intel SVML's AVX-512 high-accuracy sin breakpoint
+// tables byte-for-byte (SVML packs the fields in a different, per-precision
+// order; this table keeps one order for both precisions). Two conventions
+// are per-precision, pinned down by byte comparison against SVML:
+//   - low rounding:  float RN full precision / double RZ 24-bit
+//   - sigma cut:     bump 2^k -> 2^(k+1) when |derivative| > cut * 2^k,
+//                    cut = sqrt(2) for float / 1.5 for double
+//
+// Two structural tricks keep every rounding and comparison decidable
+// (Sollya brackets irrational values but can never prove one is exactly
+// zero, nor decide a comparison whose sides are exactly equal):
+//   - Quadrant reduction: evaluate sin/cos only in [0, pi/2) and select by
+//     quadrant. Where f or f' is exactly 0/+-1 (angle a multiple of pi/2)
+//     the reduced angle is the exact rational 0, so sin(0)/cos(0) evaluate
+//     symbolically and those entries store clean +0.0.
+//   - Cut nudge: scaling the float cut by (1 - 2^-200) resolves the exact
+//     tie |derivative| = 2^(-1/2) upward like SVML does; no other entry on
+//     either grid comes anywhere near a cut point.
 //
 // Parameters:
-//   pT        - Type descriptor with .kSize, .kDigits, .kRound
-//   pFunc     - Function to approximate (sin or cos)
-//   pFuncDriv - Derivative of pFunc (cos or -sin)
-procedure ApproxLut4_(pT, pFunc, pFuncDriv) {
+//   pT    - Type descriptor with .kSize, .kDigits, .kRound
+//   pQuad - Quadrant offset of the function to approximate:
+//           0 for sin (derivative cos), 1 for cos (derivative -sin),
+//           using cos(x) = sin(x + pi/2)
+procedure ApproxLut4_(pT, pQuad) {
   var r, i, $;
   // Table size: 512 entries for 64-bit, 256 for 32-bit
-  // More entries for double precision to maintain accuracy
   // These sizes balance table memory usage with interpolation accuracy:
   // - 256 entries = 1.4° spacing for float (sufficient for 24-bit mantissa)
   // - 512 entries = 0.7° spacing for double (needed for 53-bit mantissa)
   $.num_lut = match pT.kSize
     with 64: (2^9)
     default: (2^8);
-    
-  // Low part rounding configuration:
-  // - 64-bit: Round to 24 bits with round-to-zero (faster, sufficient for residual)
-  // - 32-bit: Use full precision with round-to-nearest
-  $.low_round = match pT.kSize 
+
+  // Low part rounding: full precision RN for float, truncated 24-bit for double
+  $.low_round = match pT.kSize
     with 64: ([|24, RZ|])
     default: ([|pT.kDigits, RN|]);
-    
-  // Scale factor to convert table index to angle in radians
-  $.scale = 2.0 * pi / $.num_lut;
-  
+
+  // Sigma cut: geometric nearest power of two for float, linear for double,
+  // nudged below the float tie point (see header)
+  $.cut = match pT.kSize
+    with 64: (3 / 2)
+    default: (sqrt(2));
+  $.cut = $.cut * (1 - 2^(-200));
+
   r = [||];
   for i from 0 to $.num_lut - 1 do {
-    // Sample angle uniformly distributed from 0 to 2π
-    $.angle = i * $.scale;
-    
-    // Compute exact function value
-    $.exact = pFunc($.angle);
-    
-    // Split into high and low parts for extended precision
-    // High part gets the main value rounded to type precision
+    // Quadrant reduction: angle = quad*pi/2 + reduced with reduced in [0, pi/2)
+    $.quadrant = floor(4 * i / $.num_lut);
+    $.reduced = 2 * pi * (i - $.quadrant * $.num_lut / 4) / $.num_lut;
+    $.sin_red = sin($.reduced);
+    $.cos_red = cos($.reduced);
+    $.sin_by_quad = [| $.sin_red, $.cos_red, -$.sin_red, -$.cos_red |];
+    $.cos_by_quad = [| $.cos_red, -$.sin_red, -$.cos_red, $.sin_red |];
+
+    // Shift by the function's quadrant offset:
+    // f(angle) = sin(angle + pQuad*pi/2), f'(angle) = cos(angle + pQuad*pi/2)
+    $.quadrant = $.quadrant + pQuad;
+    if ($.quadrant > 3) then $.quadrant = $.quadrant - 4;
+
+    // Exact function value, split into high and low parts:
+    // value ≈ high + low with high rounded to type precision
+    $.exact = $.sin_by_quad[$.quadrant];
     $.high = pT.kRound($.exact);
-    
-    // Low part stores the residual, rounded to reduced precision
-    // This allows accurate reconstruction: value ≈ high + low
     $.low = pT.kRound(round($.exact - $.high, $.low_round[0], $.low_round[1]));
-    
-    // Compute derivative for interpolation
-    $.deriv_exact = pFuncDriv($.angle);
-    
-    // Find power-of-2 scale factor closest to derivative magnitude
-    // This allows efficient storage and reconstruction
-    $.k = ceil(log2(abs($.deriv_exact)));
-    if ($.deriv_exact < 0) then $.k = -$.k;
-    
-    // Sigma is the power-of-2 scale factor
-    $.sigma = 2.0^$.k;
-    
-    // Store derivative minus sigma (typically a small value)
-    // Actual derivative = sigma + stored_deriv
-    $.deriv = pT.kRound($.deriv_exact - $.sigma);
-    
+
+    // Exact derivative for interpolation
+    $.deriv_exact = $.cos_by_quad[$.quadrant];
+
+    if ($.deriv_exact == 0) then {
+      // Derivative term vanishes (angle at pi/2 or 3pi/2 relative to f)
+      $.sigma = 0;
+      $.deriv = 0;
+    } else {
+      // Power-of-2 scale factor nearest the derivative under the cut rule;
+      // actual derivative = sigma + stored deriv
+      $.k = floor(log2(abs($.deriv_exact)));
+      $.sigma = 2.0^$.k;
+      if (abs($.deriv_exact) > $.cut * $.sigma) then $.sigma = 2 * $.sigma;
+      if ($.deriv_exact < 0) then $.sigma = -$.sigma;
+      $.deriv = pT.kRound($.deriv_exact - $.sigma);
+    };
+
     r = r @ [|$.deriv, $.sigma, $.high, $.low|];
   };
-  
+
   // Format as C array with 4 elements per table entry
   return CArrayT(pT, r, 4) @ ";";
 };
@@ -84,17 +109,17 @@ procedure ApproxLut4_(pT, pFunc, pFuncDriv) {
 Append(
   "template <typename T> inline constexpr char kSinApproxTable[] = {};",
   "template <> inline constexpr float kSinApproxTable<float>[] = ",
-  ApproxLut4_(Float32, sin(x), cos(x)),  // sin table with cos derivative
+  ApproxLut4_(Float32, 0),  // sin table with cos derivative
   "",
   "template <> inline constexpr double kSinApproxTable<double>[] = ",
-  ApproxLut4_(Float64, sin(x), cos(x)),
+  ApproxLut4_(Float64, 0),
   "",
   "template <typename T> inline constexpr char kCosApproxTable[] = {};",
   "template <> inline constexpr float kCosApproxTable<float>[] = ",
-  ApproxLut4_(Float32, cos(x), -sin(x)), // cos table with -sin derivative
+  ApproxLut4_(Float32, 1),  // cos table with -sin derivative
   "",
   "template <> inline constexpr double kCosApproxTable<double>[] = ",
-  ApproxLut4_(Float64, cos(x), -sin(x)),
+  ApproxLut4_(Float64, 1),
   ""
 );
 

npsr/trig/data/constants.h.sol

@@ -39,7 +39,9 @@ Append(
 
   "template <> inline constexpr double kPi<double, false>[] = " @
   KArray_(Float64, pi, [|RN, 24, 24, 24|], [|RN, 53|]), // no FMA: 3x24-bit + 1x53-bit
-  // The 24-bit pieces ensure n*πᵢ is exact for |n| < 2^29
+  // The 24-bit pieces ensure n*πᵢ is exact (no rounding) for the low path's
+  // largest quotient; the 53-bit tail sits at 2^-76 so its inexact product
+  // rounds ~2^-105, harmless even for tiny reduced arguments near k*π/2
   "",
 
   // Special 35-bit precision π for specific algorithms
@@ -61,18 +63,17 @@ Append(
 );
 
 // Non-FMA version of π/16 for High precision implementation
-// Special handling: components are reordered [0,2,3,1] for proper evaluation
-// Without FMA, multiplication order matters to minimize rounding errors
-vNFma = Constants(pi/16, [|RN, 27, 27|], [|RN, 29|], [|RN, 53|]);
+// Pieces are applied in this natural (descending) order; each subtraction's
+// rounding error is captured into r_lo by the reduction DAG. The 27/27/29-bit
+// heads keep n*πᵢ/16 exact for |n| <= 85445660 (= round(2^24 * 16/π), the
+// largest quotient on the high path); the 53-bit tail sits at 2^-91 so its
+// inexact product rounds ~2^-117, harmless even for tiny reduced arguments
 Append(
   "template <bool FMA> inline constexpr double kPiDiv16Prec29[] = " @
   KArray_(Float64, pi/16, [|RN, 53|], [|RN, 29|], [|RN, 53|]),
-  
-  // Non-FMA version reorders components: [0], [2], [3], [1]
-  // This ordering ensures proper evaluation without FMA:
-  // r = x - n*π₁/16 - n*π₃/16 - n*π₄/16 - n*π₂/16
+
   "template <> inline constexpr double kPiDiv16Prec29<false>[] = " @
-  CArray([|vNFma[0], vNFma[2], vNFma[3], vNFma[1]|], 4) @ ";",
+  KArray_(Float64, pi/16, [|RN, 27, 27|], [|RN, 29|], [|RN, 53|]),
   "",
 
   // Simple scalar constants