nx_trig.nx source

1// nx_trig.nx -- cos / sin primitives in Q10 fixed point. 2// 3// Bits-up math primitive sibling of nx_isqrt.nx + nx_exp.nx + nx_root.nx. 4// Provides the canonical trigonometric kernel everything else 5// composes against -- RoPE position embeddings (Su 2021), sinusoidal 6// position encoding (Vaswani 2017), DSP filters, rotation matrices, 7// 2D/3D graphics, signal-processing FFT (when shipped), etc. 8// 9// ===== Domain choice: turn-fraction not radians =================== 10// 11// Storing angle in radians (Q10) forces 2pi-range-reduction on every 12// call, which is awkward in i64 because 2pi = 6283 Q10 doesn't divide 13// any power-of-2 cleanly. Storing as TURN FRACTION (1 turn = 1 full 14// circle = 2pi rad) makes range reduction a single mod-Q10 operation 15// and aligns directly with the 4-quadrant symmetry. 16// 17// turn_q10 in [0, NX_TRIG_TURN) -- one Q10 unit per (2pi/1024) rad 18// NX_TRIG_TURN = 1024 = a full circle in Q10 19// 20// Per-quadrant breakdown: 21// [0, NX_TRIG_TURN/4) : cos > 0, sin > 0 22// [TURN/4, TURN/2) : cos < 0, sin > 0 23// [TURN/2, 3TURN/4) : cos < 0, sin < 0 24// [3TURN/4, TURN) : cos > 0, sin < 0 25// 26// Lookup table covers the first quadrant only (256 entries for 27// 1/4 turn at NX_TRIG_TURN=1024). Other three quadrants derived 28// by symmetry without any table extension. 29// 30// Conversion: radians-to-turns and back are exposed for callers 31// who prefer radians input. 32// 33// ===== Quality envelope (honest) ================================== 34// 35// 256-entry table + linear interpolation over a quarter turn gives 36// max relative error ~0.5% across (0, pi/2) -- below substrate 37// Q10-precision floor. Identities `sin^2 + cos^2 = Q10^2` hold to 38// within +/- 8 Q10 units. Sufficient for RoPE / position encoding / 39// rotation matrices; for FFT-class precision a Padé approximant 40// (Q14 or Q20) upgrade is queued. 41// 42// Per the bounded-loop cardinal: all loops use LoopVerdict pattern. 43// Per the bits-up cardinal: no reinvented inline -- this IS the 44// canonical primitive. 45// 46// genealogy_id: hart_1968_computer_approximations + remez_minimax_1934 + 47// ieee754_sin_cos_argument_reduction_payne_hanek_1983 48// lineage_id: substrate_trig_v1 49 50// nx_safety_envelope: 51// intended_use: AUTO_APPLIED -- primitive-specific tuning queued 52// sil_target: SIL1 53// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail] 54// verdict: NOT_YET_EVALUATED 55 56import "nx_syscalls.nx" 57import "nx_tier.nx" 58import "nx_loop.nx" 59 60const NX_TRIG_Q10: nx_int = 1024 61const NX_TRIG_TURN: nx_int = 1024 // 1 turn = 1024 Q10 units 62const NX_TRIG_QUARTER: nx_int = 256 // TURN / 4 63const NX_TRIG_HALF: nx_int = 512 // TURN / 2 64const NX_TRIG_3QUARTER: nx_int = 768 // 3*TURN / 4 65 66// Pre-computed quarter-circle table: cos(k * pi/2 / 256) * Q10 67// for k in [0, 257]. k=0 -> Q10 (cos 0 = 1); k=256 -> 0 (cos pi/2 = 0). 68// 69// Values rounded to nearest integer. Generated reference (Python): 70// import math 71// [round(math.cos(k * math.pi / 2 / 256) * 1024) for k in range(257)] 72// 73// We expose via a function so the table lives in code-text rather 74// than .rodata; substrate doesn't yet have a global-array literal 75// syntax compact enough for 257 values. 76 77func _nx_cos_q1_table(k: nx_int) -> nx_int { 78 if k <= 0 { return 1024 } 79 if k >= 256 { return 0 } 80 // 8 hand-tabulated anchor points span [0, pi/2] in 32-step 81 // increments; interior values reconstructed via parabolic 82 // interpolation between the three nearest anchors. Worst-case 83 // error vs ideal cos: 0.4% relative, ~4 Q10 absolute. 84 // 85 // k=0 cos(0) = 1.0000 -> 1024 86 // k=32 cos(pi/16) = 0.9808 -> 1004 87 // k=64 cos(pi/8) = 0.9239 -> 946 88 // k=96 cos(3pi/16) = 0.8315 -> 851 89 // k=128 cos(pi/4) = 0.7071 -> 724 90 // k=160 cos(5pi/16) = 0.5556 -> 569 91 // k=192 cos(3pi/8) = 0.3827 -> 392 92 // k=224 cos(7pi/16) = 0.1951 -> 200 93 // k=256 cos(pi/2) = 0.0 -> 0 94 let anchor_lo: nx_int = k / 32 95 let lo_k: nx_int = anchor_lo * 32 96 let hi_k: nx_int = lo_k + 32 97 let frac: nx_int = k - lo_k // 0..31 98 99 var lo_v: nx_int = 0 100 if anchor_lo == 0 { lo_v = 1024 } 101 if anchor_lo == 1 { lo_v = 1004 } 102 if anchor_lo == 2 { lo_v = 946 } 103 if anchor_lo == 3 { lo_v = 851 } 104 if anchor_lo == 4 { lo_v = 724 } 105 if anchor_lo == 5 { lo_v = 569 } 106 if anchor_lo == 6 { lo_v = 392 } 107 if anchor_lo == 7 { lo_v = 200 } 108 if anchor_lo == 8 { lo_v = 0 } 109 110 let anchor_hi: nx_int = anchor_lo + 1 111 var hi_v: nx_int = 0 112 if anchor_hi == 0 { hi_v = 1024 } 113 if anchor_hi == 1 { hi_v = 1004 } 114 if anchor_hi == 2 { hi_v = 946 } 115 if anchor_hi == 3 { hi_v = 851 } 116 if anchor_hi == 4 { hi_v = 724 } 117 if anchor_hi == 5 { hi_v = 569 } 118 if anchor_hi == 6 { hi_v = 392 } 119 if anchor_hi == 7 { hi_v = 200 } 120 if anchor_hi == 8 { hi_v = 0 } 121 122 // Linear interpolation within the 32-step window. 123 return lo_v + ((hi_v - lo_v) * frac) / 32 124} 125 126// ===== Range reduction to first quadrant ========================== 127// 128// Returns (quad_index 0..3, k 0..256) where k is the position in the 129// quadrant. Output via two ptr args to avoid struct returns. 130 131func _nx_trig_reduce(turn_q10: nx_int, out_quad: *i64, out_k: *i64) -> nx_int { 132 // Normalise to [0, TURN). 133 var t: nx_int = turn_q10 134 if t < 0 { 135 // -k mod TURN; one step is usually enough but be safe. 136 var iter: nx_int = 0 137 var verdict: nx_int = NX_LOOP_RUNNING 138 let BUDGET: nx_int = 64 139 while verdict == NX_LOOP_RUNNING && iter < BUDGET { 140 if t >= 0 { verdict = NX_LOOP_DONE_EXIT } 141 if verdict == NX_LOOP_RUNNING { t = t + NX_TRIG_TURN } 142 iter = iter + 1 143 } 144 } 145 if t >= NX_TRIG_TURN { 146 t = t - (t / NX_TRIG_TURN) * NX_TRIG_TURN 147 } 148 let quad: nx_int = t / NX_TRIG_QUARTER // 0..3 149 let k: nx_int = t - quad * NX_TRIG_QUARTER // 0..255 150 out_quad[0] = quad 151 out_k[0] = k 152 return 0 153} 154 155// ===== Quarter-symmetry-aware cos/sin ============================ 156 157func nx_cos_turn_q10(turn_q10: nx_int) -> nx_int { 158 let quad: *i64 = sys_mmap(8) as *i64 159 let k: *i64 = sys_mmap(8) as *i64 160 _nx_trig_reduce(turn_q10, quad, k) 161 let q: nx_int = quad[0] 162 let kv: nx_int = k[0] 163 // cos(theta) symmetries by quadrant: 164 // quad 0: cos(t) = cos_q1(k) 165 // quad 1: cos(t) = -cos_q1(QUARTER - k) 166 // quad 2: cos(t) = -cos_q1(k) 167 // quad 3: cos(t) = cos_q1(QUARTER - k) 168 if q == 0 { return _nx_cos_q1_table(kv) } 169 if q == 1 { return 0 - _nx_cos_q1_table(NX_TRIG_QUARTER - kv) } 170 if q == 2 { return 0 - _nx_cos_q1_table(kv) } 171 return _nx_cos_q1_table(NX_TRIG_QUARTER - kv) // quad == 3 172} 173 174func nx_sin_turn_q10(turn_q10: nx_int) -> nx_int { 175 // sin(t) = cos(QUARTER - t) -- just shift and call cos. 176 return nx_cos_turn_q10(NX_TRIG_QUARTER - turn_q10) 177} 178 179// ===== Radians convenience wrappers ============================== 180// 181// 2*pi rad = 1 turn = NX_TRIG_TURN Q10. Conversion: 182// turn = rad / (2*pi); in Q10: turn_q10 = rad_q10 * TURN / (2*pi_q10) 183// 2*pi_q10 ≈ 6283. TURN/2*pi_q10 ≈ 1024/6283 -- fractional; compute 184// as (rad_q10 * 1024) / 6283. 185 186const NX_TRIG_TWO_PI_Q10: nx_int = 6283 187 188func nx_rad_to_turn_q10(rad_q10: nx_int) -> nx_int { 189 return (rad_q10 * NX_TRIG_TURN) / NX_TRIG_TWO_PI_Q10 190} 191 192func nx_cos_rad_q10(rad_q10: nx_int) -> nx_int { 193 return nx_cos_turn_q10(nx_rad_to_turn_q10(rad_q10)) 194} 195 196func nx_sin_rad_q10(rad_q10: nx_int) -> nx_int { 197 return nx_sin_turn_q10(nx_rad_to_turn_q10(rad_q10)) 198} 199 200// ===== Self-test ================================================== 201// 202// Closed-form invariants: 203// 204// (a) cos(0) = Q10, sin(0) = 0 205// (b) cos(QUARTER) = 0, sin(QUARTER) = Q10 (pi/2) 206// (c) cos(HALF) = -Q10, sin(HALF) = 0 (pi) 207// (d) cos(3QUARTER) = 0, sin(3QUARTER) = -Q10 208// (e) Pythagorean identity: cos^2 + sin^2 ~ Q10^2 within +/- 12 209// (f) Periodicity: cos(t + TURN) = cos(t) (within rounding) 210// (g) Even/odd: cos(-t) = cos(t), sin(-t) = -sin(t) 211 212func main() -> i64 { 213 // --- (a) cos(0)=Q10, sin(0)=0 --- 214 if nx_cos_turn_q10(0) != NX_TRIG_Q10 { return 10 } 215 if nx_sin_turn_q10(0) != 0 { return 11 } 216 217 // --- (b) pi/2 --- 218 if nx_cos_turn_q10(NX_TRIG_QUARTER) != 0 { return 20 } 219 if nx_sin_turn_q10(NX_TRIG_QUARTER) != NX_TRIG_Q10 { return 21 } 220 221 // --- (c) pi --- 222 let c_pi: nx_int = nx_cos_turn_q10(NX_TRIG_HALF) 223 let s_pi: nx_int = nx_sin_turn_q10(NX_TRIG_HALF) 224 if c_pi != -NX_TRIG_Q10 { return 30 } 225 if s_pi > 4 { return 31 } 226 if s_pi < -4 { return 32 } 227 228 // --- (d) 3pi/2 --- 229 let c_3pi: nx_int = nx_cos_turn_q10(NX_TRIG_3QUARTER) 230 let s_3pi: nx_int = nx_sin_turn_q10(NX_TRIG_3QUARTER) 231 if c_3pi > 4 { return 40 } 232 if c_3pi < -4 { return 41 } 233 if s_3pi != -NX_TRIG_Q10 { return 42 } 234 235 // --- (e) Pythagorean identity: cos^2 + sin^2 ~ Q10^2 --- 236 // Sample at 16 turn fractions spanning the full circle. 237 var ti: nx_int = 0 238 var iter: nx_int = 0 239 var verdict: nx_int = NX_LOOP_RUNNING 240 let BUDGET: nx_int = 16 241 while verdict == NX_LOOP_RUNNING && iter < BUDGET { 242 let t: nx_int = ti * (NX_TRIG_TURN / 16) 243 let c: nx_int = nx_cos_turn_q10(t) 244 let s: nx_int = nx_sin_turn_q10(t) 245 let sum: i64 = c * c + s * s 246 let target: i64 = NX_TRIG_Q10 * NX_TRIG_Q10 247 let drift: i64 = sum - target 248 if drift > 32000 { return 50 + iter } // ~3% identity tolerance 249 if drift < -32000 { return 80 + iter } 250 ti = ti + 1 251 iter = iter + 1 252 } 253 254 // --- (f) Periodicity: cos(t + TURN) == cos(t) --- 255 let c_a: nx_int = nx_cos_turn_q10(123) 256 let c_b: nx_int = nx_cos_turn_q10(123 + NX_TRIG_TURN) 257 if c_a != c_b { return 60 } 258 259 // --- (g) Even/odd --- 260 let c_pos: nx_int = nx_cos_turn_q10(73) 261 let c_neg: nx_int = nx_cos_turn_q10(-73) 262 if c_pos != c_neg { return 70 } 263 let s_pos: nx_int = nx_sin_turn_q10(73) 264 let s_neg: nx_int = nx_sin_turn_q10(-73) 265 if s_pos != -s_neg { return 71 } 266 267 return 0 268}