Skip to content

Commit ca4f3e6

Browse files
seiko2plusseiko2plus
authored
FEAT: Python/Sollya driver that used for generates LUT/Constants (#5)
* FEAT: Python/Sollya driver that used for generates LUT/Constants Argument reduction & tables: - Payne–Hanek tables via chunked 4/π per exponent range (`reduction.h.sol`) - π splits for FMA/non-FMA and float/double (`constants.h.sol`) - Compact 4-term (value + derivatives) LUTs (`approx.h.sol`) - Precomputed sin/cos(k·π/16) for quadrant handling (`kpi16-inl.h.sol`) Sollya-driven codegen pipeline: - Tooling lives in `tools/sollya` with shared helpers - CLI: `spin sollya [-f]` (via `.spin/cmds.py`) to (re)generate headers - Project wiring: `pyproject.toml` matters for spin * 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.
1 parent b56d1b5 commit ca4f3e6

13 files changed

Lines changed: 1179 additions & 3 deletions

File tree

‎.gitignore‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
# Temporary #
2+
#############
3+
.cache
4+
.direnv
5+
6+
# Compiled source #
7+
###################
8+
*.a
9+
*.dll
10+
*.exe
11+
*.o
12+
*.o.d
13+
*.py[ocd]
14+
*.so
15+
*.mod
16+
17+
# Logs #
18+
########
19+
*.log
20+
21+
# Patches #
22+
###########
23+
*.patch
24+
*.diff

‎.spin/cmds.py‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
import pathlib
2+
import sys
3+
import click
4+
5+
curdir = pathlib.Path(__file__).parent
6+
rootdir = curdir.parent
7+
toolsdir = rootdir / "tools"
8+
sys.path.insert(0, str(toolsdir))
9+
10+
11+
@click.command(help="Generate sollya c++/python based files")
12+
@click.option("-f", "--force", is_flag=True, help="Force regenerate all files")
13+
def sollya(*, force):
14+
import sollya # type: ignore[import]
15+
16+
sollya.main(force=force)

‎npsr/trig/data/approx.h.sol‎

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
// Suppress rounding mode information messages for irrational table values
2+
suppressmessage(184, 185, 186); // suppress info no rounding, round-up, round-down
3+
suppressmessage(160); // inequality decided by faithful evaluation
4+
5+
// Generates a 4-element lookup table for fast trigonometric function approximation.
6+
//
7+
// The table format per angle is: [deriv, sigma, high, low]
8+
// where:
9+
// sigma + deriv = f'(angle) sigma = sign * 2^k nearest the derivative;
10+
// being a power of two, sigma*z is exact for any
11+
// small remainder z, confining the interpolation
12+
// rounding error to the tiny deriv*z term
13+
// high + low = f(angle) double-float split of the function value
14+
//
15+
// The values reproduce Intel SVML's AVX-512 high-accuracy sin breakpoint
16+
// tables byte-for-byte (SVML packs the fields in a different, per-precision
17+
// order; this table keeps one order for both precisions). Two conventions
18+
// are per-precision, pinned down by byte comparison against SVML:
19+
// - low rounding: float RN full precision / double RZ 24-bit
20+
// - sigma cut: bump 2^k -> 2^(k+1) when |derivative| > cut * 2^k,
21+
// cut = sqrt(2) for float / 1.5 for double
22+
//
23+
// Two structural tricks keep every rounding and comparison decidable
24+
// (Sollya brackets irrational values but can never prove one is exactly
25+
// zero, nor decide a comparison whose sides are exactly equal):
26+
// - Quadrant reduction: evaluate sin/cos only in [0, pi/2) and select by
27+
// quadrant. Where f or f' is exactly 0/+-1 (angle a multiple of pi/2)
28+
// the reduced angle is the exact rational 0, so sin(0)/cos(0) evaluate
29+
// symbolically and those entries store clean +0.0.
30+
// - Cut nudge: scaling the float cut by (1 - 2^-200) resolves the exact
31+
// tie |derivative| = 2^(-1/2) upward like SVML does; no other entry on
32+
// either grid comes anywhere near a cut point.
33+
//
34+
// Parameters:
35+
// pT - Type descriptor with .kSize, .kDigits, .kRound
36+
// pQuad - Quadrant offset of the function to approximate:
37+
// 0 for sin (derivative cos), 1 for cos (derivative -sin),
38+
// using cos(x) = sin(x + pi/2)
39+
procedure ApproxLut4_(pT, pQuad) {
40+
var r, i, $;
41+
// Table size: 512 entries for 64-bit, 256 for 32-bit
42+
// These sizes balance table memory usage with interpolation accuracy:
43+
// - 256 entries = 1.4° spacing for float (sufficient for 24-bit mantissa)
44+
// - 512 entries = 0.7° spacing for double (needed for 53-bit mantissa)
45+
$.num_lut = match pT.kSize
46+
with 64: (2^9)
47+
default: (2^8);
48+
49+
// Low part rounding: full precision RN for float, truncated 24-bit for double
50+
$.low_round = match pT.kSize
51+
with 64: ([|24, RZ|])
52+
default: ([|pT.kDigits, RN|]);
53+
54+
// Sigma cut: geometric nearest power of two for float, linear for double,
55+
// nudged below the float tie point (see header)
56+
$.cut = match pT.kSize
57+
with 64: (3 / 2)
58+
default: (sqrt(2));
59+
$.cut = $.cut * (1 - 2^(-200));
60+
61+
r = [||];
62+
for i from 0 to $.num_lut - 1 do {
63+
// Quadrant reduction: angle = quad*pi/2 + reduced with reduced in [0, pi/2)
64+
$.quadrant = floor(4 * i / $.num_lut);
65+
$.reduced = 2 * pi * (i - $.quadrant * $.num_lut / 4) / $.num_lut;
66+
$.sin_red = sin($.reduced);
67+
$.cos_red = cos($.reduced);
68+
$.sin_by_quad = [| $.sin_red, $.cos_red, -$.sin_red, -$.cos_red |];
69+
$.cos_by_quad = [| $.cos_red, -$.sin_red, -$.cos_red, $.sin_red |];
70+
71+
// Shift by the function's quadrant offset:
72+
// f(angle) = sin(angle + pQuad*pi/2), f'(angle) = cos(angle + pQuad*pi/2)
73+
$.quadrant = $.quadrant + pQuad;
74+
if ($.quadrant > 3) then $.quadrant = $.quadrant - 4;
75+
76+
// Exact function value, split into high and low parts:
77+
// value ≈ high + low with high rounded to type precision
78+
$.exact = $.sin_by_quad[$.quadrant];
79+
$.high = pT.kRound($.exact);
80+
$.low = pT.kRound(round($.exact - $.high, $.low_round[0], $.low_round[1]));
81+
82+
// Exact derivative for interpolation
83+
$.deriv_exact = $.cos_by_quad[$.quadrant];
84+
85+
if ($.deriv_exact == 0) then {
86+
// Derivative term vanishes (angle at pi/2 or 3pi/2 relative to f)
87+
$.sigma = 0;
88+
$.deriv = 0;
89+
} else {
90+
// Power-of-2 scale factor nearest the derivative under the cut rule;
91+
// actual derivative = sigma + stored deriv
92+
$.k = floor(log2(abs($.deriv_exact)));
93+
$.sigma = 2.0^$.k;
94+
if (abs($.deriv_exact) > $.cut * $.sigma) then $.sigma = 2 * $.sigma;
95+
if ($.deriv_exact < 0) then $.sigma = -$.sigma;
96+
$.deriv = pT.kRound($.deriv_exact - $.sigma);
97+
};
98+
99+
r = r @ [|$.deriv, $.sigma, $.high, $.low|];
100+
};
101+
102+
// Format as C array with 4 elements per table entry
103+
return CArrayT(pT, r, 4) @ ";";
104+
};
105+
106+
// Generate C++ header content with specialized lookup tables
107+
// for both float and double precision sine and cosine
108+
// Template declarations (empty for unsupported types)
109+
Append(
110+
"template <typename T> inline constexpr char kSinApproxTable[] = {};",
111+
"template <> inline constexpr float kSinApproxTable<float>[] = ",
112+
ApproxLut4_(Float32, 0), // sin table with cos derivative
113+
"",
114+
"template <> inline constexpr double kSinApproxTable<double>[] = ",
115+
ApproxLut4_(Float64, 0),
116+
"",
117+
"template <typename T> inline constexpr char kCosApproxTable[] = {};",
118+
"template <> inline constexpr float kCosApproxTable<float>[] = ",
119+
ApproxLut4_(Float32, 1), // cos table with -sin derivative
120+
"",
121+
"template <> inline constexpr double kCosApproxTable<double>[] = ",
122+
ApproxLut4_(Float64, 1),
123+
""
124+
);
125+
126+
WriteCPPHeader("npsr::trig::data");

‎npsr/trig/data/constants.h.sol‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
// Suppress rounding mode information messages that are expected during constant generation
2+
suppressmessage(185, 186); // suppress expected info round-up, round-down
3+
4+
// Helper procedure to format constants array for C++ output
5+
// Takes a type descriptor followed by constant values and formats them
6+
// as a C array with 4 elements per line
7+
procedure KArray_(pArgs = ...) {
8+
var pT;
9+
pT = head(pArgs);
10+
return CArrayT(pT, Constants @ tail(pArgs), 4) @ ";";
11+
};
12+
13+
// Generate C++ header with various π-related constants for Cody-Waite reduction
14+
// These constants enable accurate computation of r = x - n*π with extended precision
15+
//
16+
// The Cody-Waite method splits π into multiple parts:
17+
// r = x - n*π₁ - n*π₂ - n*π₃ ...
18+
// where each πᵢ has limited precision to ensure exact multiplication
19+
//
20+
// Different versions are provided for:
21+
// - FMA (Fused Multiply-Add) vs non-FMA architectures
22+
// - float vs double precision
23+
// - Various precision requirements
24+
Append(
25+
// Generic template declaration
26+
"template <typename T, bool FMA> inline constexpr char kPi[] = {};",
27+
28+
// Float π constants for Low precision implementation
29+
"template <> inline constexpr float kPi<float, true>[] = " @
30+
KArray_(Float32, pi, [|RN, 24, 24, 24|]), // FMA: 3x24-bit pieces (full precision each)
31+
32+
"template <> inline constexpr float kPi<float, false>[] = " @
33+
KArray_(Float32, pi, [|RD, 11, 11, 11|], [|RN, 24|]), // no FMA: 3x11-bit + 1x24-bit
34+
// The 11-bit pieces ensure n*πᵢ is exact (no rounding) for |n| < 2^13
35+
36+
// Double π constants for Low precision implementation
37+
"template <> inline constexpr double kPi<double, true>[] = " @
38+
KArray_(Float64, pi, [|RN, 53|], [|RD, 53|], [|RU, 53|]), // FMA: Different roundings for error compensation
39+
40+
"template <> inline constexpr double kPi<double, false>[] = " @
41+
KArray_(Float64, pi, [|RN, 24, 24, 24|], [|RN, 53|]), // no FMA: 3x24-bit + 1x53-bit
42+
// The 24-bit pieces ensure n*πᵢ is exact (no rounding) for the low path's
43+
// largest quotient; the 53-bit tail sits at 2^-76 so its inexact product
44+
// rounds ~2^-105, harmless even for tiny reduced arguments near k*π/2
45+
"",
46+
47+
// Special 35-bit precision π for specific algorithms
48+
"template <bool FMA> inline constexpr double kPiPrec35[] = " @
49+
KArray_(Float64, pi, [|RN, 35|], [|RD, 53|]),
50+
51+
"template <> inline constexpr double kPiPrec35<false>[] = " @
52+
KArray_(Float64, pi, [|RN, 24, 24, 24|]),
53+
"",
54+
55+
// 2π constants for angle wrapping
56+
"template <typename T> inline constexpr char kPiMul2[] = {};",
57+
"template <> inline constexpr float kPiMul2<float>[] = " @
58+
KArray_(Float32, pi*2, [|RN, 24, 24|]), // 2x24-bit pieces
59+
60+
"template <> inline constexpr double kPiMul2<double>[] = " @
61+
KArray_(Float64, pi*2, [|RN, 53, 53|]), // 2x53-bit pieces
62+
""
63+
);
64+
65+
// Non-FMA version of π/16 for High precision implementation
66+
// Pieces are applied in this natural (descending) order; each subtraction's
67+
// rounding error is captured into r_lo by the reduction DAG. The 27/27/29-bit
68+
// heads keep n*πᵢ/16 exact for |n| <= 85445660 (= round(2^24 * 16/π), the
69+
// largest quotient on the high path); the 53-bit tail sits at 2^-91 so its
70+
// inexact product rounds ~2^-117, harmless even for tiny reduced arguments
71+
Append(
72+
"template <bool FMA> inline constexpr double kPiDiv16Prec29[] = " @
73+
KArray_(Float64, pi/16, [|RN, 53|], [|RN, 29|], [|RN, 53|]),
74+
75+
"template <> inline constexpr double kPiDiv16Prec29<false>[] = " @
76+
KArray_(Float64, pi/16, [|RN, 27, 27|], [|RN, 29|], [|RN, 53|]),
77+
"",
78+
79+
// Simple scalar constants
80+
"template <typename T> inline constexpr char kInvPi = '_';",
81+
"template <> inline constexpr float kInvPi<float> = " @ single(1/pi) @ "f;",
82+
"template <> inline constexpr double kInvPi<double> = " @ double(1/pi) @ ";",
83+
"",
84+
85+
"template <typename T> inline constexpr char kHalfPi = '_';",
86+
"template <> inline constexpr float kHalfPi<float> = " @ single(pi/2) @ "f;",
87+
"template <> inline constexpr double kHalfPi<double> = " @ double(pi/2) @ ";",
88+
"",
89+
90+
"template <typename T> inline constexpr char k16DivPi = '_';",
91+
"template <> inline constexpr float k16DivPi<float> = " @ single(16/pi) @ "f;",
92+
"template <> inline constexpr double k16DivPi<double> = " @ double(16/pi) @ ";",
93+
""
94+
);
95+
// Dump();
96+
97+
WriteCPPHeader("npsr::trig::data");
98+

‎npsr/trig/data/data.h‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,11 @@
1+
// Aggregates the trig data headers. Hand-written: it only forwards includes,
2+
// so there is nothing for Sollya to generate.
3+
//
4+
// Intentionally NOT guarded with #ifndef: it pulls in the Highway target-toggled
5+
// header kpi16-inl.h, which must be re-included once per SIMD target via
6+
// hwy/foreach_target.h. An include-once guard would suppress all but the first
7+
// target pass. The include-once children (constants/approx/reduction) carry
8+
// their own guards and no-op on re-entry.
19
#include "npsr/lut-inl.h"
210
#include "npsr/trig/data/constants.h"
311
#include "npsr/trig/data/kpi16-inl.h"

‎npsr/trig/data/kpi16-inl.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ inline constexpr auto kKPi16Table = MakeLut<double>(
3535
}
3636

3737
);
38-
inline HWY_ATTR void _dummy_supress_unused_target(){}
38+
inline HWY_ATTR void _dummy_suppress_unused_target(){}
3939
} // namespace npsr::HWY_NAMESPACE::trig
4040
HWY_AFTER_NAMESPACE();
4141
#endif // NPSR_TRIG_DATA_KPI16_INL_H

‎npsr/trig/data/kpi16-inl.h.sol‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
// Generates lookup table for sin(k·π/16) and cos(k·π/16) values
2+
// Used in the high-precision trigonometric implementation for range reduction
3+
//
4+
// This table supports the algorithm where input x is reduced to:
5+
// x = n*(π/16) + r, where |r| < π/16
6+
// Then sin(x) and cos(x) are reconstructed using angle addition formulas
7+
//
8+
// Parameters:
9+
// pT - Type descriptor (Float64 in this case)
10+
// pFunc - Function to evaluate (sin or cos)
11+
// pBy - Divisor for π (16 in this case, giving π/16 intervals)
12+
procedure PiDivTable_(pT, pFunc, pBy) {
13+
var r, i, pi_by;
14+
pi_by = pi / pBy;
15+
r = [||];
16+
17+
// Generate function values at k*π/16 for k = 0, 1, ..., 15
18+
for i from 0 to pBy - 1 do {
19+
r = r :. pT.kRound(pFunc(i * pi_by));
20+
};
21+
22+
// Format as C array with 4 elements per line
23+
return CArrayT(pT, r, 4);
24+
};
25+
26+
// Generates packed low-precision parts of sin and cos values
27+
// This packing scheme saves memory by storing two 32-bit values in one 64-bit word
28+
//
29+
// The packing works as follows for double (64-bit):
30+
// - sin_low occupies bits [31:0] (lower 32 bits)
31+
// - cos_low occupies bits [63:32] (upper 32 bits)
32+
//
33+
// This is why in the C++ code:
34+
// - cos_lo can be used directly (it's already in the upper bits)
35+
// - sin_lo needs to be extracted with a 32-bit left shift
36+
procedure PiDivPackLowTable_(pT, pFunc0, pFunc1, pBy) {
37+
var r, i, digits, $;
38+
$.pi_by = pi / pBy;
39+
r = [||];
40+
41+
// First, compute the low precision parts (residuals after high precision)
42+
for i from 0 to pBy - 1 do {
43+
$.hi0 = pT.kRound(pFunc0(i * $.pi_by)); // High precision sin
44+
$.hi1 = pT.kRound(pFunc1(i * $.pi_by)); // High precision cos
45+
// Low precision parts: exact value minus high precision part
46+
$.hi0_low = pT.kRound(pFunc0(i * $.pi_by) - $.hi0); // sin_low
47+
$.hi1_low = pT.kRound(pFunc1(i * $.pi_by) - $.hi1); // cos_low
48+
r = r @ [|$.hi0_low, $.hi1_low|];
49+
};
50+
51+
// Convert to binary representation for bit manipulation
52+
digits = ToDigits(pT, r);
53+
$.half_size = pT.kSize / 2; // 32 for double
54+
$.lower_bits = 2^$.half_size; // Mask for lower 32 bits
55+
56+
r = [||];
57+
// Pack pairs of values into single 64-bit words
58+
for i from 0 to length(digits) - 1 by 2 do {
59+
$.hi0 = digits[i]; // sin_low bits
60+
$.hi1 = digits[i + 1]; // cos_low bits
61+
$.pack = mod(RightShift($.hi0, $.half_size), $.lower_bits);
62+
$.pack = $.pack + $.hi1 - mod($.hi1, $.lower_bits);
63+
r = r :. $.pack;
64+
};
65+
66+
// Convert back from binary representation
67+
r = FromDigits(pT, r);
68+
return CArrayT(pT, r, 4);
69+
};
70+
71+
Append(
72+
"inline constexpr auto kKPi16Table = MakeLut<double>(",
73+
"// High parts of sin(k·π/16) where k = 0, 1, ..., 15",
74+
PiDivTable_(Float64, sin(x), 16) @ ",",
75+
"// High parts of cos(k·π/16) where k = 0, 1, ..., 15",
76+
PiDivTable_(Float64, cos(x), 16) @ ",",
77+
"// Lower parts of sin(k·π/16) and cos(k·π/16) packed together",
78+
"// Format: bits [63:32] = cos_low, bits [31:0] = sin_low",
79+
"// This packing saves 16×8 = 128 bytes of memory",
80+
PiDivPackLowTable_(Float64, sin(x), cos(x), 16),
81+
"",
82+
");"
83+
);
84+
85+
WriteHighwayHeader("npsr::HWY_NAMESPACE::trig");
86+

‎npsr/trig/data/reduction.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -262,7 +262,7 @@ template <> inline constexpr uint32_t kLargeReductionTable<float>[] = {
262262
2795223067u, 1817359143u, 2289174591u,
263263
1295478838u, 3634718287u, 283381887u,
264264
2590957677u, 2974469278u, 566763775u,
265-
};;
265+
};
266266

267267
template <> inline constexpr uint64_t kLargeReductionTable<double>[] = {
268268
0ull, 0ull, 0ull,
@@ -2313,7 +2313,7 @@ template <> inline constexpr uint64_t kLargeReductionTable<double>[] = {
23132313
2862820329136421693ull, 520719405400910443ull, 18113390975381211222ull,
23142314
5725640658272843386ull, 1041438810801820887ull, 17780037877052870828ull,
23152315
11451281316545686772ull, 2082877621603641775ull, 17113331680396190040ull,
2316-
};;
2316+
};
23172317

23182318
} // namespace npsr::trig::data
23192319

0 commit comments

Comments
 (0)