nx_i128.nx source

1// nx_i128.nx -- minimum-viable signed 128-bit integer primitive. 2// 3// Numeric tier N2 in the substrate's tier ladder (see docs/ 4// NUMERIC_TIER_LADDER.md). Exists because i64 (N1) hits overflow 5// in cases that occur in the wild: e.g. HLL cubic Lagrange 6// interpolation at lg_k=10 with milli-units X table produces 7// intermediate products up to ~1.7e21, exceeding i64 max ~9.2e18. 8// 9// The substrate's discipline: when i64 isn't enough, ESCALATE to 10// i128 (or higher). Never silently downgrade precision. 11// 12// Representation: { hi: i64, lo: i64 } where lo holds the low 64 13// bits (treated as unsigned via masking) and hi holds the upper 64 14// bits (signed). Two's-complement throughout. Sign of the i128 15// = sign of hi. 16// 17// Operations: enough to support Lagrange interpolation 18// (mul i64*i64 -> i128, shl by k, div by i64 -> i64). More 19// operations added as additional callers need them. 20// 21// Numerical Recipes ref: §1.1 (basic arithmetic) + §5.6 (Lagrange). 22// Hacker's Delight ch.8 (signed/unsigned multiword arithmetic). 23// genealogy_id: knuth_taocp_vol2 + warren_hackers_delight 24// lineage_id: integer_arithmetic 25 26// nx_safety_envelope: 27// intended_use: AUTO_APPLIED -- primitive-specific tuning queued 28// sil_target: SIL1 29// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail] 30// verdict: NOT_YET_EVALUATED 31 32import "syscalls.nx" 33 34struct I128 { 35 hi: i64, 36 lo: i64, 37} 38 39const NX_I128_LO_MASK: i64 = 0xFFFFFFFF 40const NX_I128_HALF_BITS: i64 = 32 41 42// === construction ==================================================== 43 44func nx_i128_alloc() -> *I128 { 45 let raw: *u8 = sys_mmap(16) 46 let z: *I128 = raw as *I128 47 z.hi = 0 48 z.lo = 0 49 return z 50} 51 52func nx_i128_set_i64(z: *I128, v: i64) -> i64 { 53 z.lo = v 54 if v < 0 { z.hi = -1 } 55 if v >= 0 { z.hi = 0 } 56 return 0 57} 58 59// === sign / negation ================================================= 60 61func nx_i128_is_neg(z: *I128) -> i64 { 62 if z.hi < 0 { return 1 } 63 return 0 64} 65 66// Two's-complement negation: invert all bits, add 1. 67// NishiLang lacks the `~` unary; use `x ^ -1` (since -1 has all bits set). 68func nx_i128_neg(z: *I128) -> i64 { 69 let inv_lo: i64 = z.lo ^ -1 70 let inv_hi: i64 = z.hi ^ -1 71 let new_lo: i64 = inv_lo + 1 72 var new_hi: i64 = inv_hi 73 // Carry from lo: when original z.lo == 0, inv_lo == -1, +1 wraps to 0 74 // and we propagate +1 to hi. 75 if z.lo == 0 { new_hi = new_hi + 1 } 76 z.lo = new_lo 77 z.hi = new_hi 78 return 0 79} 80 81// === add / sub ======================================================= 82 83func nx_i128_add(z: *I128, b: *I128) -> i64 { 84 let a_lo_u: i64 = z.lo & NX_I128_LO_MASK 85 let a_lo_h: i64 = (z.lo >> NX_I128_HALF_BITS) & NX_I128_LO_MASK 86 let b_lo_u: i64 = b.lo & NX_I128_LO_MASK 87 let b_lo_h: i64 = (b.lo >> NX_I128_HALF_BITS) & NX_I128_LO_MASK 88 let sum_u: i64 = a_lo_u + b_lo_u 89 let carry_u: i64 = sum_u >> NX_I128_HALF_BITS 90 let sum_h: i64 = a_lo_h + b_lo_h + carry_u 91 let carry_h: i64 = sum_h >> NX_I128_HALF_BITS 92 z.lo = (sum_u & NX_I128_LO_MASK) | ((sum_h & NX_I128_LO_MASK) << NX_I128_HALF_BITS) 93 z.hi = z.hi + b.hi + carry_h 94 return 0 95} 96 97// === unsigned 64*64 -> 128 multiply ================================== 98// 99// Splits each i64 into two 32-bit unsigned halves and computes four 100// 32x32 -> 64 partial products that all fit in i64 (each up to ~4.3e9 101// * 4.3e9 = 1.8e19... wait that overflows. But (2^32-1)^2 = 2^64 - 102// 2^33 + 1 ~ 1.8e19 — DOES overflow i64. Need to be careful.) 103// 104// Actually (2^32 - 1)^2 = 0xFFFFFFFE_00000001 which is > 2^63. So 105// even a 32x32 multiply can overflow signed i64. We treat it as 106// unsigned, knowing the bits are correct even if the i64 interp is 107// "negative". Subsequent shifts/masks recover the right halves. 108 109func nx_i128_mul_u64(a: i64, b: i64, out: *I128) -> i64 { 110 let a_lo: i64 = a & NX_I128_LO_MASK 111 let a_hi: i64 = (a >> NX_I128_HALF_BITS) & NX_I128_LO_MASK 112 let b_lo: i64 = b & NX_I128_LO_MASK 113 let b_hi: i64 = (b >> NX_I128_HALF_BITS) & NX_I128_LO_MASK 114 115 let ll: i64 = a_lo * b_lo 116 let lh: i64 = a_lo * b_hi 117 let hl: i64 = a_hi * b_lo 118 let hh: i64 = a_hi * b_hi 119 120 // Assemble: prod = ll + (lh + hl) << 32 + hh << 64 121 // ll <= (2^32-1)^2 ~ 1.8e19 -- treated as unsigned 64-bit. 122 let ll_lo: i64 = ll & NX_I128_LO_MASK 123 let ll_hi: i64 = (ll >> NX_I128_HALF_BITS) & NX_I128_LO_MASK 124 let mid: i64 = lh + hl + ll_hi // up to 3 * (2^32-1)^2 / 2^32 ~ 1.3e10, fits 125 let mid_lo: i64 = mid & NX_I128_LO_MASK 126 let mid_hi: i64 = (mid >> NX_I128_HALF_BITS) & NX_I128_LO_MASK 127 out.lo = ll_lo | (mid_lo << NX_I128_HALF_BITS) 128 out.hi = hh + mid_hi 129 return 0 130} 131 132// === signed 64*64 -> 128 multiply ==================================== 133 134func nx_i128_mul_i64(a: i64, b: i64, out: *I128) -> i64 { 135 var neg: i64 = 0 136 var ua: i64 = a 137 var ub: i64 = b 138 if a < 0 { ua = -a; neg = neg ^ 1 } 139 if b < 0 { ub = -b; neg = neg ^ 1 } 140 nx_i128_mul_u64(ua, ub, out) 141 if neg == 1 { nx_i128_neg(out) } 142 return 0 143} 144 145// === shl ============================================================= 146// 147// Left shift in place. n must be in [0, 127]. 148 149func nx_i128_shl(z: *I128, n: i64) -> i64 { 150 if n == 0 { return 0 } 151 if n >= 64 { 152 z.hi = z.lo << (n - 64) 153 z.lo = 0 154 return 0 155 } 156 let lo_top: i64 = (z.lo >> (64 - n)) & ((1 << n) - 1) 157 z.hi = (z.hi << n) | lo_top 158 z.lo = z.lo << n 159 return 0 160} 161 162// === unsigned i128 / i64 = i64 (with quotient assumed to fit i64) ==== 163// 164// Long division by 32-bit halves. Knuth Algorithm D simplified for 165// the 32-bit-divisor case; we extend to 64-bit divisor by splitting 166// the divisor as well. Sufficient for the HLL Lagrange case where 167// numerator is i128 but result is known i64. 168// 169// This is the "schoolbook" 128/64 -> 64 routine. Two-step: 170// 1. Compute hi_dividend = z.hi (treated unsigned) 171// 2. Long-divide by d, accumulating quotient bits. 172// 256 iterations of shift-and-subtract. 173 174func nx_i128_udiv_u64(z: *I128, d: i64) -> i64 { 175 if d == 0 { return 0 } // caller responsibility 176 // Special case: hi == 0 -> simple 64-bit divide. 177 if z.hi == 0 { 178 if z.lo >= 0 { return z.lo / d } 179 // z.lo treated unsigned but i64 negative bit-pattern means 180 // value >= 2^63. Split: lo = (high_part * 2 + low_bit) ... 181 // For our HLL use case the dividend in i128 always has hi > 0 182 // when overflow occurred, so z.lo < 0 with z.hi == 0 won't 183 // occur in practice. Conservative path: shift-and-subtract. 184 } 185 // General shift-and-subtract using bit-by-bit long division. 186 var rem_hi: i64 = 0 187 var rem_lo: i64 = 0 188 var q_hi: i64 = 0 189 var q_lo: i64 = 0 190 var bit: i64 = 127 191 while bit >= 0 { 192 // Shift remainder left 1. 193 let new_rem_hi: i64 = (rem_hi << 1) | ((rem_lo >> 63) & 1) 194 let new_rem_lo: i64 = rem_lo << 1 195 rem_hi = new_rem_hi 196 rem_lo = new_rem_lo 197 // Bring in next dividend bit. 198 var src_bit: i64 = 0 199 if bit >= 64 { 200 src_bit = (z.hi >> (bit - 64)) & 1 201 } 202 if bit < 64 { 203 src_bit = (z.lo >> bit) & 1 204 } 205 rem_lo = rem_lo | src_bit 206 // Compare remainder to d (unsigned). rem_hi==0 means rem<2^64, 207 // can compare rem_lo directly with d when d>0. 208 // For our use case where final quotient fits i64, rem_hi will 209 // be 0 by the time we're computing meaningful bits. 210 var ge: i64 = 0 211 if rem_hi > 0 { ge = 1 } 212 if rem_hi == 0 { 213 if rem_lo < 0 { ge = 1 } // rem_lo as unsigned >= 2^63 > d (d > 0, d < 2^63) 214 if rem_lo >= 0 { 215 if rem_lo >= d { ge = 1 } 216 } 217 } 218 if ge == 1 { 219 // rem -= d 220 let sub_lo: i64 = rem_lo - d 221 var borrow: i64 = 0 222 // Detect borrow on the unsigned subtract: borrow if (rem_lo as unsigned) < d. 223 if rem_lo >= 0 { 224 if rem_lo < d { borrow = 1 } 225 } 226 // If rem_lo < 0 (as i64), it represents unsigned >= 2^63 >= d, no borrow. 227 rem_lo = sub_lo 228 rem_hi = rem_hi - borrow 229 // Set quotient bit. 230 if bit >= 64 { 231 q_hi = q_hi | (1 << (bit - 64)) 232 } 233 if bit < 64 { 234 q_lo = q_lo | (1 << bit) 235 } 236 } 237 bit = bit - 1 238 } 239 return q_lo // we assume quotient fits i64 240} 241 242// Signed i128 / i64 -> i64. 243func nx_i128_div_i64(z: *I128, d: i64) -> i64 { 244 if d == 0 { return 0 } 245 var neg: i64 = 0 246 if nx_i128_is_neg(z) == 1 { 247 nx_i128_neg(z) 248 neg = neg ^ 1 249 } 250 var ud: i64 = d 251 if d < 0 { ud = -d; neg = neg ^ 1 } 252 let q: i64 = nx_i128_udiv_u64(z, ud) 253 if neg == 1 { return -q } 254 return q 255} 256 257// === convenience: muldiv (a*b/c via i128) =========================== 258// 259// The HLL Lagrange workhorse. Computes (a * b) / c with i128 260// intermediate so we don't lose precision when a*b overflows i64. 261// Returns i64 quotient (caller ensures result fits). 262 263func nx_muldiv_i64(a: i64, b: i64, c: i64) -> i64 { 264 let prod: *I128 = nx_i128_alloc() 265 nx_i128_mul_i64(a, b, prod) 266 let q: i64 = nx_i128_div_i64(prod, c) 267 return q 268} 269 270// === convenience: mulshldiv ((a*b) << shift) / c via i128 ========== 271// 272// Variant for Q-format ratios where we want (a*b * SCALE) / c. 273// Computes a*b in i128, shifts left by shift bits (still i128), 274// divides by c, returns i64. 275 276func nx_mulshl_div_i64(a: i64, b: i64, shift: i64, c: i64) -> i64 { 277 let prod: *I128 = nx_i128_alloc() 278 nx_i128_mul_i64(a, b, prod) 279 nx_i128_shl(prod, shift) 280 let q: i64 = nx_i128_div_i64(prod, c) 281 return q 282}