nx_f64_oracle_sov.nx source
↩ module page · 373 lines · 12591 B
1// nx_f64_oracle_sov.nx -- SOVEREIGN f64 oracle: pure-Nishi independent recompute
2// of IEEE 754 binary64 add/sub/mul/div/sqrt expected values (operator law: no py,
3// no sh -- the Python vector generator becomes a retired debt for ME0).
4//
5// INDEPENDENCE: this is a structurally different path from nx_f64. The impl
6// rounds through per-op guard/sticky shortcuts (sticky-jam GRS); the oracle
7// carries the EXACT result in a two-word 120-bit significand (hi:lo, 60 bits
8// per word) plus a sticky flag, finds the true msb, and rounds from the full
9// value. Anchor proof: _f64_soak_gate_authored phase 1 reproduces ALL 1217
10// hardware-IEEE-anchored KAT expected values; phase 2 soaks random vectors
11// impl-vs-oracle (sovereign xorshift PRNG).
12//
13// Word convention: value = (hi*2^60 + lo) * 2^(e112 - 1135); a normal 53-bit
14// significand sig placed at hi=sig, lo=0 has its msb at bit 112 and biased
15// exponent e112. orc_pack derives eo = e112 + (msb - 112), denormalizes below
16// eo=1 into sticky, rounds to nearest even from the exact tail.
17//
18// license_tier: ORIGINAL
19
20import "nx_syscalls.nx"
21import "nx_tier.nx"
22const ORC_MAGIC_2047: i64 = 2047
23const ORC_MAGIC_1075: i64 = 1075
24const ORC_MAGIC_1107: i64 = 1107
25
26const ORC_M60: i64 = 0x0FFFFFFFFFFFFFFF // 2^60 - 1
27const ORC_NAN: i64 = 0x7FF8000000000000
28const ORC_INF: i64 = 0x7FF0000000000000
29const ORC_IMP: i64 = 0x0010000000000000 // 1 << 52
30const ORC_M52: i64 = 0x000FFFFFFFFFFFFF
31
32func orc_sign(raw: i64) -> i64 { return (raw >> 63) & 1 }
33func orc_ef(raw: i64) -> i64 { return (raw >> 52) & 0x7FF }
34func orc_mf(raw: i64) -> i64 { return raw & ORC_M52 }
35
36// class: 0 zero, 1 normal, 2 subnormal, 3 inf, 4 nan
37func orc_cls(raw: i64) -> i64 {
38 let e: i64 = orc_ef(raw)
39 let m: i64 = orc_mf(raw)
40 if e == 0 { if m == 0 { return 0 } return 2 }
41 if e == ORC_MAGIC_2047 { if m == 0 { return 3 } return 4 }
42 return 1
43}
44
45// significand normalized to [2^52, 2^53); returns sig, writes effective biased
46// exponent through *eout (subnormals shift left, exponent goes <= 0).
47func orc_sig_norm(raw: i64, eout: *i64) -> i64 {
48 var sig: i64 = orc_mf(raw)
49 var ef: i64 = orc_ef(raw)
50 if ef == 0 {
51 ef = 1
52 while sig < ORC_IMP { sig = sig << 1; ef = ef - 1 }
53 } else {
54 sig = sig | ORC_IMP
55 }
56 eout[0] = ef
57 return sig
58}
59
60// bit i of (hi:lo), i in [0, 122]
61func orc_bit(hi: i64, lo: i64, i: i64) -> i64 {
62 if i < 60 { return (lo >> i) & 1 }
63 return (hi >> (i - 60)) & 1
64}
65
66// any set bit strictly below position i?
67func orc_below(hi: i64, lo: i64, i: i64) -> i64 {
68 if i <= 0 { return 0 }
69 if i <= 60 {
70 if (lo & ((1 << i) - 1)) != 0 { return 1 }
71 return 0
72 }
73 if lo != 0 { return 1 }
74 let k: i64 = i - 60
75 if k >= 63 { if hi != 0 { return 1 } return 0 }
76 if (hi & ((1 << k) - 1)) != 0 { return 1 }
77 return 0
78}
79
80// 53 bits of (hi:lo) starting at bit position s (s >= 0; bits above msb are 0).
81func orc_extract53(hi: i64, lo: i64, s: i64) -> i64 {
82 if s >= 60 { return hi >> (s - 60) }
83 var m: i64 = lo >> s
84 let need_hi: i64 = s - 7 // number of hi bits that land in the window
85 if need_hi > 0 {
86 m = m | ((hi & ((1 << need_hi) - 1)) << (60 - s))
87 }
88 // need_hi <= 0 implies the msb sits inside lo (hi == 0): nothing to merge.
89 return m & ((1 << 53) - 1)
90}
91
92// Round the exact value (hi:lo, sticky) * 2^(e112-1135) to nearest-even f64.
93func orc_pack(sign: i64, e112: i64, hi: i64, lo: i64, sticky_in: i64) -> i64 {
94 var sticky: i64 = sticky_in
95 if hi == 0 {
96 if lo == 0 {
97 return sign << 63 // exact zero (sticky can only round-to-zero here)
98 }
99 }
100 // msb position
101 var p: i64 = 0
102 if hi != 0 {
103 var t: i64 = hi
104 var k: i64 = 60
105 while t > 1 { t = t >> 1; k = k + 1 }
106 p = k
107 } else {
108 var t2: i64 = lo
109 var k2: i64 = 0
110 while t2 > 1 { t2 = t2 >> 1; k2 = k2 + 1 }
111 p = k2
112 }
113 var eo: i64 = e112 + (p - 112)
114 var shift: i64 = p - 52
115 // denormalize: keep eo >= 1 by extracting from a higher shift
116 if eo < 1 {
117 shift = shift + (1 - eo)
118 eo = 1
119 }
120 var mant: i64 = 0
121 var guard: i64 = 0
122 if shift <= 0 {
123 // small exact value (deep cancellation): left shift, no rounding bits
124 mant = lo << (0 - shift)
125 if hi != 0 { mant = mant | (hi << (60 - shift)) }
126 } else {
127 if shift > 122 {
128 if orc_below(hi, lo, 123) != 0 { sticky = 1 }
129 mant = 0
130 guard = 0
131 } else {
132 mant = orc_extract53(hi, lo, shift)
133 guard = orc_bit(hi, lo, shift - 1)
134 if orc_below(hi, lo, shift - 1) != 0 { sticky = 1 }
135 }
136 }
137 // round to nearest even
138 var up: i64 = 0
139 if guard == 1 {
140 if sticky == 1 { up = 1 }
141 if up == 0 { if (mant & 1) == 1 { up = 1 } }
142 }
143 if up == 1 {
144 mant = mant + 1
145 if mant >= (1 << 53) { mant = mant >> 1; eo = eo + 1 }
146 }
147 if eo >= ORC_MAGIC_2047 { if mant >= ORC_IMP { return (sign << 63) | ORC_INF } }
148 if mant < ORC_IMP { return (sign << 63) | mant }
149 return (sign << 63) | (eo << 52) | (mant & ORC_M52)
150}
151
152// ===== addition / subtraction (subtract = add of negated b) ==========
153
154func orc_add(a: i64, b: i64) -> i64 {
155 let ca: i64 = orc_cls(a)
156 let cb: i64 = orc_cls(b)
157 if ca == 4 { return ORC_NAN }
158 if cb == 4 { return ORC_NAN }
159 let sa: i64 = orc_sign(a)
160 let sb: i64 = orc_sign(b)
161 if ca == 3 {
162 if cb == 3 { if sa == sb { return a } return ORC_NAN }
163 return a
164 }
165 if cb == 3 { return b }
166 if ca == 0 {
167 if cb == 0 { if sa == sb { return a } return 0 }
168 return b
169 }
170 if cb == 0 { return a }
171
172 let ea_p: *i64 = sys_mmap(16) as *i64
173 let eb_p: *i64 = sys_mmap(16) as *i64
174 var siga: i64 = orc_sig_norm(a, ea_p)
175 var sigb: i64 = orc_sig_norm(b, eb_p)
176 var ea: i64 = ea_p[0]
177 var eb: i64 = eb_p[0]
178 var s_a: i64 = sa
179 var s_b: i64 = sb
180 if eb > ea {
181 let t1: i64 = siga; siga = sigb; sigb = t1
182 let t2: i64 = ea; ea = eb; eb = t2
183 let t3: i64 = s_a; s_a = s_b; s_b = t3
184 }
185 let diff: i64 = ea - eb
186 // big at (hi=siga, lo=0); small shifted right by diff across the 120-bit frame
187 var sh: i64 = 0
188 var sl: i64 = 0
189 var st: i64 = 0
190 if diff <= 60 {
191 sh = sigb >> diff
192 if diff == 0 {
193 sl = 0
194 } else {
195 sl = (sigb & ((1 << diff) - 1)) << (60 - diff)
196 }
197 } else {
198 if diff <= 113 {
199 sh = 0
200 sl = sigb >> (diff - 60)
201 let dlow: i64 = diff - 60
202 if dlow < 63 {
203 if (sigb & ((1 << dlow) - 1)) != 0 { st = 1 }
204 }
205 } else {
206 sh = 0
207 sl = 0
208 st = 1
209 }
210 }
211
212 var rh: i64 = 0
213 var rl: i64 = 0
214 var rs: i64 = s_a
215 if s_a == s_b {
216 rl = sl
217 rh = siga + sh
218 // (big lo word is 0, so no lo carry)
219 } else {
220 // (siga : 0) - (sh : sl), then if sticky bits were lost from small the
221 // true small is a touch LARGER: decrement one more lsb and keep sticky
222 // (the value lies in (D-1, D); representing as (D-1)+sticky rounds right).
223 var borrow: i64 = 0
224 if sl > 0 {
225 rl = (ORC_M60 + 1) - sl
226 borrow = 1
227 } else {
228 rl = 0
229 }
230 rh = siga - sh - borrow
231 if st == 1 {
232 if rl == 0 { rl = ORC_M60; rh = rh - 1 } else { rl = rl - 1 }
233 }
234 if rh < 0 {
235 // |b| > |a|: only possible at diff == 0 (then sl = 0, st = 0): flip.
236 rh = sh - siga
237 rl = 0
238 rs = s_b
239 }
240 if rh == 0 { if rl == 0 { if st == 0 { return 0 } } }
241 }
242 return orc_pack(rs, ea, rh, rl, st)
243}
244
245func orc_sub(a: i64, b: i64) -> i64 {
246 return orc_add(a, b ^ (1 << 63))
247}
248
249// ===== multiplication =================================================
250
251func orc_mul(a: i64, b: i64) -> i64 {
252 let ca: i64 = orc_cls(a)
253 let cb: i64 = orc_cls(b)
254 let so: i64 = orc_sign(a) ^ orc_sign(b)
255 if ca == 4 { return ORC_NAN }
256 if cb == 4 { return ORC_NAN }
257 if ca == 3 { if cb == 0 { return ORC_NAN } return (so << 63) | ORC_INF }
258 if cb == 3 { if ca == 0 { return ORC_NAN } return (so << 63) | ORC_INF }
259 if ca == 0 { return so << 63 }
260 if cb == 0 { return so << 63 }
261
262 let ea_p: *i64 = sys_mmap(16) as *i64
263 let eb_p: *i64 = sys_mmap(16) as *i64
264 let siga: i64 = orc_sig_norm(a, ea_p)
265 let sigb: i64 = orc_sig_norm(b, eb_p)
266
267 // exact 106-bit product via 20/33-bit split (different split than the impl)
268 let a0: i64 = siga & ((1 << 33) - 1)
269 let a1: i64 = siga >> 33 // 20 bits
270 let b0: i64 = sigb & ((1 << 33) - 1)
271 let b1: i64 = sigb >> 33
272 let hh: i64 = a1 * b1 // <= 2^40, at bit 66
273 let hl: i64 = a1 * b0 + a0 * b1 // < 2^55, at bit 33
274 let ll: i64 = a0 * b0 // < 2^66 -- CAREFUL: fits? 33+33=66 > 63!
275 // a0,b0 < 2^33 -> product < 2^66 overflows i64. Split ll differently:
276 // use a0 = a0h*2^16 + a0l to keep partials < 2^50.
277 let a0h: i64 = a0 >> 16
278 let a0l: i64 = a0 & 0xFFFF
279 let p0: i64 = a0h * b0 // < 2^50, at bit 16
280 let p1: i64 = a0l * b0 // < 2^49, at bit 0
281 // assemble into lo (bits 0..59) and hi (bits 60..)
282 // contributions: p1@0, p0@16, hl@33, hh@66
283 let lo1: i64 = p1 & ORC_M60
284 let c1: i64 = p1 >> 60
285 let p0_lo: i64 = (p0 & ((1 << 44) - 1)) << 16
286 let p0_hi: i64 = p0 >> 44
287 let hl_lo: i64 = (hl & ((1 << 27) - 1)) << 33
288 let hl_hi: i64 = hl >> 27
289 var lo: i64 = lo1 + p0_lo + hl_lo
290 var carry: i64 = lo >> 60
291 lo = lo & ORC_M60
292 let hi: i64 = (hh << 6) + p0_hi + hl_hi + c1 + carry
293
294 // value = P * 2^(ea+eb-2150); P = hi*2^60+lo -> e112 = ea+eb-1015
295 return orc_pack(so, ea_p[0] + eb_p[0] - 1015, hi, lo, 0)
296}
297
298// ===== division =======================================================
299
300func orc_div(a: i64, b: i64) -> i64 {
301 let ca: i64 = orc_cls(a)
302 let cb: i64 = orc_cls(b)
303 let so: i64 = orc_sign(a) ^ orc_sign(b)
304 if ca == 4 { return ORC_NAN }
305 if cb == 4 { return ORC_NAN }
306 if cb == 0 { if ca == 0 { return ORC_NAN } return (so << 63) | ORC_INF }
307 if cb == 3 { if ca == 3 { return ORC_NAN } return so << 63 }
308 if ca == 0 { return so << 63 }
309 if ca == 3 { return (so << 63) | ORC_INF }
310
311 let ea_p: *i64 = sys_mmap(16) as *i64
312 let eb_p: *i64 = sys_mmap(16) as *i64
313 let siga: i64 = orc_sig_norm(a, ea_p)
314 let sigb: i64 = orc_sig_norm(b, eb_p)
315
316 // q = floor(siga * 2^60 / sigb) bit-serial; r exact
317 var q: i64 = 0
318 var r: i64 = siga
319 if r >= sigb { q = 1; r = r - sigb }
320 var i: i64 = 0
321 while i < 60 {
322 r = r << 1
323 q = q << 1
324 if r >= sigb { q = q + 1; r = r - sigb }
325 i = i + 1
326 }
327 var st: i64 = 0
328 if r != 0 { st = 1 }
329 // value = q * 2^(ea-eb-60+...) -> e112 = ea - eb + 1075
330 let hi: i64 = q >> 60
331 let lo: i64 = q & ORC_M60
332 return orc_pack(so, ea_p[0] - eb_p[0] + ORC_MAGIC_1075, hi, lo, st)
333}
334
335// ===== square root ====================================================
336
337func orc_sqrt(a: i64) -> i64 {
338 let ca: i64 = orc_cls(a)
339 if ca == 4 { return ORC_NAN }
340 if ca == 0 { return a }
341 if orc_sign(a) == 1 { return ORC_NAN }
342 if ca == 3 { return a }
343
344 let e_p: *i64 = sys_mmap(16) as *i64
345 var sig: i64 = orc_sig_norm(a, e_p)
346 var e_unb: i64 = e_p[0] - 1023
347 if (e_unb & 1) != 0 { sig = sig << 1; e_unb = e_unb - 1 }
348 let f_half: i64 = (e_unb - 52) >> 1
349
350 // root = floor(sqrt(sig * 2^56)) in [2^54, 2^55): 28 pairs of sig (56 bits)
351 // then 28 zero pairs
352 var rem: i64 = 0
353 var root: i64 = 0
354 var idx: i64 = 27
355 while idx >= 0 {
356 let pair: i64 = (sig >> (idx * 2)) & 3
357 rem = (rem << 2) | pair
358 let trial: i64 = (root << 2) | 1
359 if rem >= trial { rem = rem - trial; root = (root << 1) | 1 } else { root = root << 1 }
360 idx = idx - 1
361 }
362 var j: i64 = 0
363 while j < 28 {
364 rem = rem << 2
365 let trial2: i64 = (root << 2) | 1
366 if rem >= trial2 { rem = rem - trial2; root = (root << 1) | 1 } else { root = root << 1 }
367 j = j + 1
368 }
369 var st: i64 = 0
370 if rem != 0 { st = 1 }
371 // value = root * 2^(f_half - 28) -> e112 = f_half + 1107
372 return orc_pack(0, f_half + ORC_MAGIC_1107, 0, root, st)
373}