nx_numeric.nx source

1// numeric.nx -- Kahan-discipline numerical primitives. 2// 3// Per W. Kahan ("How Java's Floating-Point Hurts Everyone Everywhere" 4// 1998 + "MxMulEps" matrix-mul precision paper), a numerically sound 5// language must expose substrate-level building blocks for 6// catastrophic-cancellation-resistant arithmetic. This file ships 7// the canonical three: TwoSum, TwoProduct, and Kahan compensated 8// summation. 9// 10// Compilation note: this module uses `f64` types and arithmetic 11// operators directly, which the NishiLang sovereign sibling 12// (nxc.elf) accepts but the C bootstrap (nxc2.exe) does NOT lex. 13// That is the correct division of labor: C bootstrap is one 14// Wheeler-comparator anchor for the integer subset; f64-using 15// modules Wheeler against a different peer (rustc, gcc-riscv64 16// compiling an equivalent C program, or self-host nxc.elf when 17// it ships). NishiLang is not limited by what its bootstrap 18// can lex; the bootstrap is limited by NishiLang's needs. 19// 20// Calling convention: state structs caller-allocated via 21// sys_mmap, mutated through pointers; scalar helpers return 22// single f64 values. No struct-by-value returns, no &local 23// out-params (pre-existing language constraints). 24// 25// Bibliography: 26// Dekker (1971) "A Floating-Point Technique for Extending 27// the Available Precision" -- TwoSum, 28// TwoProduct, the SPLIT_FACTOR trick 29// Møller (1965) pre-Dekker accumulating-correction summation 30// Kahan (1965) compensated summation algorithm 31// Kahan & Ivory (1995) "Roundoff Degrades an Idealized Cantilever" 32// Ogita-Rump-Oishi (2005) "Accurate Sum and Dot Product" 33// -- compensated dot product 34 35// nx_safety_envelope: 36// intended_use: AUTO_APPLIED -- primitive-specific tuning queued 37// sil_target: SIL1 38// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail] 39// verdict: NOT_YET_EVALUATED 40 41import "nx_syscalls.nx" 42 43// === Kahan running-sum state ====================================== 44// 45// 16 bytes: { s: f64, err: f64 }. s is the rounded running sum; 46// err is the pending compensation that captures bits IEEE rounding 47// has dropped so far. Reusable via nx_kahan_reset. 48 49struct KahanState { 50 s: f64, 51 err: f64, 52} 53 54func nx_kahan_alloc() -> *KahanState { 55 let raw: *u8 = sys_mmap(16) 56 let st: *KahanState = raw as *KahanState 57 st.s = 0.0 58 st.err = 0.0 59 return st 60} 61 62func nx_kahan_reset(st: *KahanState) -> i64 { 63 st.s = 0.0 64 st.err = 0.0 65 return 0 66} 67 68func nx_kahan_finish(st: *KahanState) -> f64 { 69 return st.s 70} 71 72// === TwoSum (Knuth 1969 / Møller 1965) ============================ 73// 74// Given two f64 inputs a, b, compute the IEEE-rounded sum s = a + b 75// and the EXACT rounding error err such that a + b == s + err in 76// real arithmetic. Their pair carries a > 53-bit-precision 77// representation of the true sum. 78// 79// The "fast" 6-op TwoSum that doesn't require |a| >= |b|. Dekker's 80// 3-op variant exists but needs the magnitude precondition; we 81// prefer fewer branches over fewer ops. 82// 83// API split into nx_two_sum (returns s) and nx_two_sum_err 84// (returns err) so each is single-return. Callers needing both 85// call both -- the IEEE rounding is deterministic, so s computed 86// twice gives the same bits. 87 88func nx_two_sum(a: f64, b: f64) -> f64 { 89 return a + b 90} 91 92func nx_two_sum_err(a: f64, b: f64) -> f64 { 93 let s: f64 = a + b 94 let bp: f64 = s - a 95 let ap: f64 = s - bp 96 let db: f64 = b - bp 97 let da: f64 = a - ap 98 return da + db 99} 100 101// === Dekker split (helper for TwoProduct) ========================= 102// 103// SPLIT_FACTOR = 2^27 + 1 = 134217729. Multiplying an f64 by this 104// rounds to a_hi exactly because the trailing 26 bits get carried 105// into a_hi during round-to-nearest-even; a - a_hi yields a_lo, 106// the dropped low half. 107 108const NX_DEKKER_SPLIT: f64 = 134217729.0 109 110func nx_dekker_hi(a: f64) -> f64 { 111 let c: f64 = NX_DEKKER_SPLIT * a 112 return c - (c - a) 113} 114 115func nx_dekker_lo(a: f64) -> f64 { 116 let c: f64 = NX_DEKKER_SPLIT * a 117 let hi: f64 = c - (c - a) 118 return a - hi 119} 120 121// === TwoProduct (Veltkamp / Dekker) =============================== 122// 123// Compute the IEEE-rounded product p = a * b and the exact 124// rounding error err such that a * b == p + err in real arithmetic. 125// 17 f64 ops without FMA. When OP_FMADD codegen lands the err 126// computation collapses to fma(a, b, -p) (2 ops); this version 127// uses Dekker's split for portability. 128 129func nx_two_product(a: f64, b: f64) -> f64 { 130 return a * b 131} 132 133func nx_two_product_err(a: f64, b: f64) -> f64 { 134 let p: f64 = a * b 135 let a_hi: f64 = nx_dekker_hi(a) 136 let a_lo: f64 = nx_dekker_lo(a) 137 let b_hi: f64 = nx_dekker_hi(b) 138 let b_lo: f64 = nx_dekker_lo(b) 139 let t1: f64 = a_hi * b_hi - p 140 let t2: f64 = a_hi * b_lo 141 let t3: f64 = a_lo * b_hi 142 let t4: f64 = a_lo * b_lo 143 return t1 + t2 + t3 + t4 144} 145 146// === Kahan Compensated Summation (1965) =========================== 147// 148// Add x to the running state st with O(1) extra precision: tracks 149// the per-step rounding loss in st.err and reintroduces it on the 150// next step. Backward error bound: 151// |sum - true| <= 2 * eps * sum_of_|x_i| 152// independent of n -- vs naive summation's O(n*eps) bound. 153// 154// Caller: allocate via nx_kahan_alloc, feed via nx_kahan_add, 155// retrieve via nx_kahan_finish. 156 157func nx_kahan_add(st: *KahanState, x: f64) -> i64 { 158 let y: f64 = x - st.err 159 let t: f64 = st.s + y 160 let ts: f64 = t - st.s 161 let new_err: f64 = ts - y 162 st.s = t 163 st.err = new_err 164 return 0 165} 166 167// === Compensated Dot Product (Ogita-Rump-Oishi 2005) ============== 168// 169// dot(x, y) = sum_i x_i * y_i with a parallel-error track that 170// absorbs both per-product and per-sum rounding. Backward error 171// bound matches a 2K-bit accumulator at every step. Cost: ~25 172// f64 ops per element vs naive's 2. Kahan's MxMulEps paper makes 173// this non-negotiable for trustworthy matrix multiplication. 174 175func nx_compensated_dot(x: *f64, y: *f64, n: i64) -> f64 { 176 let sum_p: *KahanState = nx_kahan_alloc() 177 let sum_e: *KahanState = nx_kahan_alloc() 178 var i: i64 = 0 179 while i < n { 180 let xi: f64 = x[i] 181 let yi: f64 = y[i] 182 let p: f64 = xi * yi 183 let e: f64 = nx_two_product_err(xi, yi) 184 nx_kahan_add(sum_p, p) 185 nx_kahan_add(sum_e, e) 186 i = i + 1 187 } 188 let final_e: f64 = sum_e.s + sum_e.err 189 nx_kahan_add(sum_p, final_e) 190 return nx_kahan_finish(sum_p) 191}