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 35import "syscalls.nx" 36 37// === Kahan running-sum state ====================================== 38// 39// 16 bytes: { s: f64, err: f64 }. s is the rounded running sum; 40// err is the pending compensation that captures bits IEEE rounding 41// has dropped so far. Reusable via nx_kahan_reset. 42 43struct KahanState { 44 s: f64, 45 err: f64, 46} 47 48func nx_kahan_alloc() -> *KahanState { 49 let raw: *u8 = sys_mmap(16) 50 let st: *KahanState = raw as *KahanState 51 st.s = 0.0 52 st.err = 0.0 53 return st 54} 55 56func nx_kahan_reset(st: *KahanState) -> i64 { 57 st.s = 0.0 58 st.err = 0.0 59 return 0 60} 61 62func nx_kahan_finish(st: *KahanState) -> f64 { 63 return st.s 64} 65 66// === TwoSum (Knuth 1969 / Møller 1965) ============================ 67// 68// Given two f64 inputs a, b, compute the IEEE-rounded sum s = a + b 69// and the EXACT rounding error err such that a + b == s + err in 70// real arithmetic. Their pair carries a > 53-bit-precision 71// representation of the true sum. 72// 73// The "fast" 6-op TwoSum that doesn't require |a| >= |b|. Dekker's 74// 3-op variant exists but needs the magnitude precondition; we 75// prefer fewer branches over fewer ops. 76// 77// API split into nx_two_sum (returns s) and nx_two_sum_err 78// (returns err) so each is single-return. Callers needing both 79// call both -- the IEEE rounding is deterministic, so s computed 80// twice gives the same bits. 81 82func nx_two_sum(a: f64, b: f64) -> f64 { 83 return a + b 84} 85 86func nx_two_sum_err(a: f64, b: f64) -> f64 { 87 let s: f64 = a + b 88 let bp: f64 = s - a 89 let ap: f64 = s - bp 90 let db: f64 = b - bp 91 let da: f64 = a - ap 92 return da + db 93} 94 95// === Dekker split (helper for TwoProduct) ========================= 96// 97// SPLIT_FACTOR = 2^27 + 1 = 134217729. Multiplying an f64 by this 98// rounds to a_hi exactly because the trailing 26 bits get carried 99// into a_hi during round-to-nearest-even; a - a_hi yields a_lo, 100// the dropped low half. 101 102const NX_DEKKER_SPLIT: f64 = 134217729.0 103 104func nx_dekker_hi(a: f64) -> f64 { 105 let c: f64 = NX_DEKKER_SPLIT * a 106 return c - (c - a) 107} 108 109func nx_dekker_lo(a: f64) -> f64 { 110 let c: f64 = NX_DEKKER_SPLIT * a 111 let hi: f64 = c - (c - a) 112 return a - hi 113} 114 115// === TwoProduct (Veltkamp / Dekker) =============================== 116// 117// Compute the IEEE-rounded product p = a * b and the exact 118// rounding error err such that a * b == p + err in real arithmetic. 119// 17 f64 ops without FMA. When OP_FMADD codegen lands the err 120// computation collapses to fma(a, b, -p) (2 ops); this version 121// uses Dekker's split for portability. 122 123func nx_two_product(a: f64, b: f64) -> f64 { 124 return a * b 125} 126 127func nx_two_product_err(a: f64, b: f64) -> f64 { 128 let p: f64 = a * b 129 let a_hi: f64 = nx_dekker_hi(a) 130 let a_lo: f64 = nx_dekker_lo(a) 131 let b_hi: f64 = nx_dekker_hi(b) 132 let b_lo: f64 = nx_dekker_lo(b) 133 let t1: f64 = a_hi * b_hi - p 134 let t2: f64 = a_hi * b_lo 135 let t3: f64 = a_lo * b_hi 136 let t4: f64 = a_lo * b_lo 137 return t1 + t2 + t3 + t4 138} 139 140// === Kahan Compensated Summation (1965) =========================== 141// 142// Add x to the running state st with O(1) extra precision: tracks 143// the per-step rounding loss in st.err and reintroduces it on the 144// next step. Backward error bound: 145// |sum - true| <= 2 * eps * sum_of_|x_i| 146// independent of n -- vs naive summation's O(n*eps) bound. 147// 148// Caller: allocate via nx_kahan_alloc, feed via nx_kahan_add, 149// retrieve via nx_kahan_finish. 150 151func nx_kahan_add(st: *KahanState, x: f64) -> i64 { 152 let y: f64 = x - st.err 153 let t: f64 = st.s + y 154 let ts: f64 = t - st.s 155 let new_err: f64 = ts - y 156 st.s = t 157 st.err = new_err 158 return 0 159} 160 161// === Compensated Dot Product (Ogita-Rump-Oishi 2005) ============== 162// 163// dot(x, y) = sum_i x_i * y_i with a parallel-error track that 164// absorbs both per-product and per-sum rounding. Backward error 165// bound matches a 2K-bit accumulator at every step. Cost: ~25 166// f64 ops per element vs naive's 2. Kahan's MxMulEps paper makes 167// this non-negotiable for trustworthy matrix multiplication. 168 169func nx_compensated_dot(x: *f64, y: *f64, n: i64) -> f64 { 170 let sum_p: *KahanState = nx_kahan_alloc() 171 let sum_e: *KahanState = nx_kahan_alloc() 172 var i: i64 = 0 173 while i < n { 174 let xi: f64 = x[i] 175 let yi: f64 = y[i] 176 let p: f64 = xi * yi 177 let e: f64 = nx_two_product_err(xi, yi) 178 nx_kahan_add(sum_p, p) 179 nx_kahan_add(sum_e, e) 180 i = i + 1 181 } 182 let final_e: f64 = sum_e.s + sum_e.err 183 nx_kahan_add(sum_p, final_e) 184 return nx_kahan_finish(sum_p) 185}