nx_bigfloat120.nx source
↩ module page · 360 lines · 11974 B
1// nx_bigfloat120.nx -- ME3 rung 1: SOVEREIGN fixed-120-bit binary float, pure Nishi.
2//
3// The in-house high-precision oracle that retires the last python/mpmath debt in the
4// math-engine arc: 120 significand bits = 67 guard bits beyond binary64, far past what
5// correctly-rounded f64 transcendental reference values need. Rung 2 (arbitrary limbs
6// on bigint.nx) extends this; the API is chosen so that swap is additive.
7//
8// Representation (struct-free, mmap blocks of 3 i64):
9// b[0] = e -- unbiased exponent; value = (hi*2^60 + lo) * 2^(e - 119)
10// b[1] = hi -- significand bits 119..60, normalized: hi in [2^59, 2^60)
11// b[2] = lo -- significand bits 59..0
12// Values are POSITIVE-ONLY (plus a zero sentinel e = BF_ZE); callers track signs at
13// the boundaries (series with all-positive terms by construction: exp uses the
14// reciprocal trick for negative r, atanh is odd so |s| drives the series). This
15// keeps the core branch-light and audit-friendly.
16//
17// Internal sub-ulp bits are truncated (no sticky): each op is exact to 2^-119
18// relative; the transcendental pipelines run < 100 ops, so accumulated error stays
19// below 2^-110 -- 57 bits of margin over the half-ulp-of-f64 decision boundary.
20//
21// Division lives in nx_bigfloat120_div.nx (same-file mul+div codegen quirk law).
22// license_tier: ORIGINAL
23
24import "nx_syscalls.nx"
25import "nx_tier.nx"
26const BF_MAGIC_2047: i64 = 2047
27
28const BF_M60: i64 = 0x0FFFFFFFFFFFFFFF // 2^60 - 1
29const BF_M30: i64 = 0x3FFFFFFF // 2^30 - 1
30const BF_HI_MIN: i64 = 0x0800000000000000 // 2^59
31const BF_HI_TOP: i64 = 0x1000000000000000 // 2^60
32const BF_ZE: i64 = 0 - 999999 // zero sentinel exponent
33
34func bf_new() -> *i64 {
35 let b: *i64 = sys_mmap(24) as *i64
36 b[0] = BF_ZE; b[1] = 0; b[2] = 0
37 return b
38}
39
40func bf_is_zero(b: *i64) -> i64 {
41 if b[0] == BF_ZE { return 1 }
42 return 0
43}
44
45func bf_copy(dst: *i64, src: *i64) -> i64 {
46 dst[0] = src[0]; dst[1] = src[1]; dst[2] = src[2]
47 return 0
48}
49
50// set from a small positive integer (1..2^53)
51func bf_set_int(b: *i64, v: i64) -> i64 {
52 if v <= 0 { b[0] = BF_ZE; b[1] = 0; b[2] = 0; return 0 }
53 var m: i64 = v
54 var k: i64 = 0
55 while m > 1 { m = m >> 1; k = k + 1 } // k = msb index
56 // place msb at bit 119: value = v * 2^(119-k) * 2^(e-119) with e = k
57 var hi: i64 = 0
58 var lo: i64 = 0
59 if k >= 60 {
60 hi = v >> (k - 59)
61 lo = (v & ((1 << (k - 59)) - 1)) << (119 - k)
62 lo = lo & BF_M60
63 } else {
64 hi = (v << (59 - k))
65 lo = 0
66 }
67 b[0] = k; b[1] = hi; b[2] = lo
68 return 0
69}
70
71// normalize a raw (e, hi, lo) where hi:lo may be de-normalized but nonzero,
72// hi < 2^62 (small overflow allowed); fixes hi into [2^59, 2^60).
73func bf_norm(b: *i64) -> i64 {
74 var e: i64 = b[0]
75 var hi: i64 = b[1]
76 var lo: i64 = b[2]
77 if hi == 0 { if lo == 0 { b[0] = BF_ZE; return 0 } }
78 while hi >= BF_HI_TOP {
79 lo = (lo >> 1) | ((hi & 1) << 59)
80 hi = hi >> 1
81 e = e + 1
82 }
83 while hi < BF_HI_MIN {
84 hi = (hi << 1) | ((lo >> 59) & 1)
85 lo = (lo << 1) & BF_M60
86 e = e - 1
87 }
88 b[0] = e; b[1] = hi; b[2] = lo
89 return 0
90}
91
92// compare |a| vs |b|: -1 / 0 / 1
93func bf_cmp(a: *i64, b: *i64) -> i64 {
94 let za: i64 = bf_is_zero(a)
95 let zb: i64 = bf_is_zero(b)
96 if za == 1 { if zb == 1 { return 0 } return 0 - 1 }
97 if zb == 1 { return 1 }
98 if a[0] < b[0] { return 0 - 1 }
99 if a[0] > b[0] { return 1 }
100 if a[1] < b[1] { return 0 - 1 }
101 if a[1] > b[1] { return 1 }
102 if a[2] < b[2] { return 0 - 1 }
103 if a[2] > b[2] { return 1 }
104 return 0
105}
106
107// two-word right shift by k (k >= 0), truncating
108func bf_shr2(hilo: *i64, k: i64) -> i64 {
109 var hi: i64 = hilo[0]
110 var lo: i64 = hilo[1]
111 if k >= 120 { hilo[0] = 0; hilo[1] = 0; return 0 }
112 if k >= 60 {
113 lo = hi >> (k - 60)
114 hi = 0
115 } else {
116 if k > 0 {
117 lo = (lo >> k) | ((hi & ((1 << k) - 1)) << (60 - k))
118 hi = hi >> k
119 }
120 }
121 hilo[0] = hi; hilo[1] = lo
122 return 0
123}
124
125// out = a + b (both positive)
126func bf_add(out: *i64, a: *i64, b: *i64) -> i64 {
127 if bf_is_zero(a) == 1 { return bf_copy(out, b) }
128 if bf_is_zero(b) == 1 { return bf_copy(out, a) }
129 // order into scalars (no pointer swap: i64 locals only)
130 var ea: i64 = a[0]; var ha: i64 = a[1]; var la: i64 = a[2]
131 var eb: i64 = b[0]; var hb: i64 = b[1]; var lb: i64 = b[2]
132 if eb > ea {
133 let t0: i64 = ea; ea = eb; eb = t0
134 let t1: i64 = ha; ha = hb; hb = t1
135 let t2: i64 = la; la = lb; lb = t2
136 }
137 let diff: i64 = ea - eb
138 let sm: *i64 = sys_mmap(16) as *i64
139 sm[0] = hb; sm[1] = lb
140 bf_shr2(sm, diff)
141 var lo: i64 = la + sm[1]
142 var carry: i64 = lo >> 60
143 lo = lo & BF_M60
144 let hi: i64 = ha + sm[0] + carry
145 out[0] = ea; out[1] = hi; out[2] = lo
146 return bf_norm(out)
147}
148
149// out = a - b; REQUIRES a >= b (caller compares first). Exact-cancel -> zero.
150func bf_sub(out: *i64, a: *i64, b: *i64) -> i64 {
151 if bf_is_zero(b) == 1 { return bf_copy(out, a) }
152 let diff: i64 = a[0] - b[0]
153 let sm: *i64 = sys_mmap(16) as *i64
154 sm[0] = b[1]; sm[1] = b[2]
155 bf_shr2(sm, diff)
156 var lo: i64 = a[2] - sm[1]
157 var borrow: i64 = 0
158 if lo < 0 { lo = lo + BF_HI_TOP; borrow = 1 }
159 let hi: i64 = a[1] - sm[0] - borrow
160 out[0] = a[0]; out[1] = hi; out[2] = lo
161 if hi == 0 { if lo == 0 { out[0] = BF_ZE; return 0 } }
162 return bf_norm(out)
163}
164
165// out = a * b (both positive, normalized). 240-bit product via 30-bit limbs,
166// top 120 kept.
167func bf_mul(out: *i64, a: *i64, b: *i64) -> i64 {
168 if bf_is_zero(a) == 1 { return bf_copy(out, a) }
169 if bf_is_zero(b) == 1 { return bf_copy(out, b) }
170 let al: *i64 = sys_mmap(32) as *i64
171 let bl: *i64 = sys_mmap(32) as *i64
172 al[0] = a[2] & BF_M30; al[1] = a[2] >> 30; al[2] = a[1] & BF_M30; al[3] = a[1] >> 30
173 bl[0] = b[2] & BF_M30; bl[1] = b[2] >> 30; bl[2] = b[1] & BF_M30; bl[3] = b[1] >> 30
174 let col: *i64 = sys_mmap(64) as *i64 // 8 column accumulators
175 var i: i64 = 0
176 while i < 4 {
177 var j: i64 = 0
178 while j < 4 {
179 let p: i64 = al[i] * bl[j]
180 col[i + j] = col[i + j] + (p & BF_M30)
181 col[i + j + 1] = col[i + j + 1] + (p >> 30)
182 j = j + 1
183 }
184 i = i + 1
185 }
186 // carry-propagate 30-bit columns
187 var c: i64 = 0
188 var k: i64 = 0
189 while k < 8 {
190 col[k] = col[k] + c
191 c = col[k] >> 30
192 col[k] = col[k] & BF_M30
193 k = k + 1
194 }
195 // product bits 0..239 in col[0..7]; msb at 239 or 238.
196 var hi: i64 = (col[7] << 30) | col[6]
197 var lo: i64 = (col[5] << 30) | col[4]
198 var e: i64 = a[0] + b[0] + 1
199 if hi < BF_HI_MIN {
200 hi = (hi << 1) | ((lo >> 59) & 1)
201 lo = (lo << 1) & BF_M60
202 // pull the next product bit (bit 119 of the discarded tail = col[3] bit 29)
203 lo = lo | ((col[3] >> 29) & 1)
204 e = e - 1
205 }
206 out[0] = e; out[1] = hi; out[2] = lo
207 return 0
208}
209
210// out = a * n for small positive integer n (n < 2^33: limb*n < 2^63; covers
211// k <= 1100 in k*ln2 and k < 2^21 in k*pio2 trig reduction)
212func bf_mul_small(out: *i64, a: *i64, n: i64) -> i64 {
213 if bf_is_zero(a) == 1 { return bf_copy(out, a) }
214 if n == 0 { out[0] = BF_ZE; out[1] = 0; out[2] = 0; return 0 }
215 let l0: i64 = (a[2] & BF_M30) * n
216 let l1: i64 = (a[2] >> 30) * n
217 let l2: i64 = (a[1] & BF_M30) * n
218 let l3: i64 = (a[1] >> 30) * n
219 var c0: i64 = l0 & BF_M30
220 var c1: i64 = l1 + (l0 >> 30)
221 var c2: i64 = l2 + (c1 >> 30)
222 c1 = c1 & BF_M30
223 var c3: i64 = l3 + (c2 >> 30)
224 c2 = c2 & BF_M30
225 let c4: i64 = c3 >> 30 // overflow limb, up to ~12 bits
226 c3 = c3 & BF_M30
227 var hi: i64 = (c3 << 30) | c2
228 var lo: i64 = (c1 << 30) | c0
229 var e: i64 = a[0]
230 if c4 > 0 {
231 // shift right by the bit-length of c4 so its msb lands at bit 119
232 var s: i64 = 0
233 var t: i64 = c4
234 while t > 0 { t = t >> 1; s = s + 1 }
235 lo = (lo >> s) | ((hi & ((1 << s) - 1)) << (60 - s))
236 hi = (hi >> s) | (c4 << (60 - s))
237 e = e + s
238 }
239 out[0] = e
240 out[1] = hi
241 out[2] = lo
242 return bf_norm(out)
243}
244
245// out = a / n for small positive integer n (n < 2^20)
246func bf_div_small(out: *i64, a: *i64, n: i64) -> i64 {
247 if bf_is_zero(a) == 1 { return bf_copy(out, a) }
248 var r: i64 = 0
249 let q3: i64 = (a[1] >> 30)
250 var t: i64 = q3
251 let d3: i64 = t / n
252 r = t - d3 * n
253 t = (r << 30) | (a[1] & BF_M30)
254 let d2: i64 = t / n
255 r = t - d2 * n
256 t = (r << 30) | (a[2] >> 30)
257 let d1: i64 = t / n
258 r = t - d1 * n
259 t = (r << 30) | (a[2] & BF_M30)
260 let d0: i64 = t / n
261 r = t - d0 * n
262 // one extra 30-bit digit of quotient from the remainder (keeps precision
263 // after the left renormalization that division by n>1 forces)
264 t = r << 30
265 let dx: i64 = t / n
266 var hi: i64 = (d3 << 30) | d2
267 var lo: i64 = (d1 << 30) | d0
268 out[0] = a[0]
269 out[1] = hi
270 out[2] = lo
271 if hi == 0 { if lo == 0 { out[0] = BF_ZE; return 0 } }
272 bf_norm(out)
273 // merge the extension digit into the (shifted-left) low end: after norm the
274 // significand moved left by s bits; the top s bits of dx belong in lo.
275 let s: i64 = a[0] - out[0]
276 if s > 0 {
277 if s <= 30 {
278 out[2] = out[2] | ((dx >> (30 - s)) & ((1 << s) - 1))
279 }
280 }
281 return 0
282}
283
284// set from positive f64 bit pattern (raw must be finite, nonzero, sign ignored)
285func bf_set_f64(b: *i64, raw: i64) -> i64 {
286 var sig: i64 = raw & 0x000FFFFFFFFFFFFF
287 var ef: i64 = (raw >> 52) & 0x7FF
288 if ef == 0 {
289 ef = 1
290 while sig < 0x0010000000000000 { sig = sig << 1; ef = ef - 1 }
291 } else {
292 sig = sig | 0x0010000000000000
293 }
294 b[0] = ef - 1023
295 b[1] = sig << 7
296 b[2] = 0
297 return 0
298}
299
300// round to f64 bit pattern with round-to-nearest-even; sign supplied; eadj added
301// to the exponent (the 2^k rescale from range reduction). Handles overflow ->
302// inf, gradual underflow -> subnormal/zero.
303func bf_to_f64(b: *i64, sign: i64, eadj: i64) -> i64 {
304 if bf_is_zero(b) == 1 { return sign << 63 }
305 let e: i64 = b[0] + eadj
306 var mant: i64 = b[1] >> 7
307 var guard: i64 = (b[1] >> 6) & 1
308 var sticky: i64 = 0
309 if (b[1] & 63) != 0 { sticky = 1 }
310 if b[2] != 0 { sticky = 1 }
311 var eo: i64 = e + 1023
312 if eo < 1 {
313 let k: i64 = 1 - eo
314 if k > 54 { return sign << 63 }
315 var i: i64 = 0
316 while i < k {
317 if guard == 1 { sticky = 1 }
318 guard = mant & 1
319 mant = mant >> 1
320 i = i + 1
321 }
322 eo = 1
323 }
324 var up: i64 = 0
325 if guard == 1 {
326 if sticky == 1 { up = 1 }
327 if up == 0 { if (mant & 1) == 1 { up = 1 } }
328 }
329 if up == 1 {
330 mant = mant + 1
331 if mant >= 0x0020000000000000 { mant = mant >> 1; eo = eo + 1 }
332 }
333 if eo >= BF_MAGIC_2047 { if mant >= 0x0010000000000000 { return (sign << 63) | 0x7FF0000000000000 } }
334 if mant < 0x0010000000000000 { return (sign << 63) | mant }
335 return (sign << 63) | (eo << 52) | (mant & 0x000FFFFFFFFFFFFF)
336}
337
338// nearest integer of a positive bf (value < 2^40); half rounds away from zero
339func bf_to_int_nearest(b: *i64) -> i64 {
340 if bf_is_zero(b) == 1 { return 0 }
341 let e: i64 = b[0]
342 if e < 0 - 1 { return 0 } // < 1/4 -> 0 (e = -1 means [0.5, 1))
343 var ip: i64 = 0
344 var half: i64 = 0
345 if e >= 0 {
346 // integer part = top (e+1) bits of the significand
347 if e <= 59 {
348 ip = b[1] >> (59 - e)
349 half = (b[1] >> (58 - e)) & 1
350 } else {
351 // values that large never occur in range reduction (k <= ~1100)
352 ip = 0
353 }
354 } else {
355 // e == -1: value in [0.5, 1): ip = 0, half = top bit (always 1)
356 ip = 0
357 half = 1
358 }
359 return ip + half
360}