code wiki / _hdl_build / _f64_pow_gate_authored.nx
_f64_pow_gate_authored.nx source
↩ module page · 200 lines · 8104 B
1// _f64_pow_gate_authored.nx -- ME1 pow gate: dd kernel (_pe_f64pow) vs the LIVE
2// sovereign bigfloat pow oracle (nx_bigfloat120_pow). Pure Nishi end to end.
3// 2D banded pairs across the hard regions: x near 1 with amplifying y (the
4// double-double stress case), the overflow seam (z ~ 1024) and the subnormal/
5// underflow seam (z ~ -1074), integer-y parity paths with negative x,
6// subnormal x, generic mid-range, and the IEEE special matrix.
7// ULP via the signed ordered map. GREEN iff max ULP <= 1.
8// Markers: PWBAD x= y= exp= got= ulp= / POW-ULP ... / POW-GATE verdict=
9
10import "nx_syscalls.nx"
11import "nx_f64.nx"
12import "nx_f64_div.nx"
13import "nx_f64_cvt.nx"
14import "nx_bigfloat120.nx"
15import "nx_bigfloat120_div.nx"
16import "nx_bigfloat120_exp.nx"
17import "nx_bigfloat120_ln.nx"
18import "nx_bigfloat120_pow.nx"
19import "_pm_pow_spec.nx"
20import "_pe_f64pow.nx"
21
22func pwg_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
23func pwg_putn(v: i64) -> i64 { let bb: *u8 = sys_mmap(28); var m: i64 = v; if m < 0 { m = 0 - m; sys_write(1, "-" as *u8, 1) }; let t: *u8 = sys_mmap(28); var k: i64 = 0; if m == 0 { t[0] = 48; k = 1 }; while m > 0 { t[k] = 48 + (m % 10); m = m / 10; k = k + 1 }; var i: i64 = 0; while i < k { bb[i] = t[k-1-i]; i = i + 1 }; sys_write(1, bb, k); return 0 }
24func pwg_puthex(v: i64) -> i64 { let bb: *u8 = sys_mmap(20); var i: i64 = 0; while i < 16 { let nib: i64 = (v >> ((15 - i) * 4)) & 15; if nib < 10 { bb[i] = 48 + nib } else { bb[i] = 55 + nib } i = i + 1 } sys_write(1, bb, 16); return 0 }
25
26func pwg_rng(state: *i64) -> i64 {
27 var s: i64 = state[0]
28 s = s ^ ((s >> 12) & 0x000FFFFFFFFFFFFF)
29 s = s ^ (s << 25)
30 s = s ^ ((s >> 27) & 0x0000001FFFFFFFFF)
31 state[0] = s
32 return s * 2685821657736338717
33}
34
35func pwg_ord(v: i64) -> i64 {
36 if v >= 0 { return v }
37 return (1 << 63) - v
38}
39
40func pwg_ulp(got: i64, want: i64) -> i64 {
41 if got == want { return 0 }
42 var u: i64 = pwg_ord(got) - pwg_ord(want)
43 if u < 0 { u = 0 - u }
44 return u
45}
46
47func main() -> i64 {
48 let xs: *i64 = sys_mmap(8 * 1280) as *i64
49 let ys: *i64 = sys_mmap(8 * 1280) as *i64
50 var n: i64 = 0
51 let one: i64 = 0x3FF0000000000000
52 let two: i64 = 0x4000000000000000
53 let half: i64 = 0x3FE0000000000000
54 let three: i64 = 0x4008000000000000
55 // ---- special matrix pairs (kernel must MATCH the oracle bit-for-bit) ----
56 xs[n] = 0x7FF8000000000000; ys[n] = 0; n = n + 1
57 xs[n] = one; ys[n] = 0x7FF8000000000000; n = n + 1
58 xs[n] = 0x7FF8000000000000; ys[n] = one; n = n + 1
59 xs[n] = 0xBFF0000000000000; ys[n] = 0x7FF0000000000000; n = n + 1
60 xs[n] = two; ys[n] = 0x7FF0000000000000; n = n + 1
61 xs[n] = two; ys[n] = 0xFFF0000000000000; n = n + 1
62 xs[n] = half; ys[n] = 0xFFF0000000000000; n = n + 1
63 xs[n] = 0x7FF0000000000000; ys[n] = three; n = n + 1
64 xs[n] = 0xFFF0000000000000; ys[n] = three; n = n + 1
65 xs[n] = 0xFFF0000000000000; ys[n] = 0xC008000000000000; n = n + 1
66 xs[n] = 0; ys[n] = three; n = n + 1
67 xs[n] = 1 << 63; ys[n] = three; n = n + 1
68 xs[n] = 1 << 63; ys[n] = 0xC008000000000000; n = n + 1
69 xs[n] = 0xC000000000000000; ys[n] = three; n = n + 1
70 xs[n] = 0xC000000000000000; ys[n] = two; n = n + 1
71 xs[n] = 0xC000000000000000; ys[n] = half; n = n + 1
72 xs[n] = 0xBFF0000000000000; ys[n] = three; n = n + 1
73 xs[n] = 0x4010000000000000; ys[n] = half; n = n + 1
74 xs[n] = two; ys[n] = 0x4090000000000000; n = n + 1 // 2^1024 = inf seam
75 xs[n] = two; ys[n] = 0xC090C80000000000; n = n + 1 // 2^-1074 floor
76 xs[n] = two; ys[n] = 0xC090CC0000000000; n = n + 1 // 2^-1075 tie->0
77 xs[n] = half; ys[n] = 0x3FF0000000000001; n = n + 1 // y just past 1
78 let st: *i64 = sys_mmap(16) as *i64
79 st[0] = 17320508075688772935
80 var j: i64 = 0
81 // ---- band A: generic x in (2^-3, 2^3), y in +-(0..128) ----
82 while j < 300 {
83 var rx: i64 = pwg_rng(st)
84 let efx: i64 = 1020 + (rx & 7)
85 xs[n] = (efx << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
86 var ry: i64 = pwg_rng(st)
87 let efy: i64 = 1023 + (ry & 7)
88 var yv: i64 = (efy << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
89 if (pwg_rng(st) & 1) == 1 { yv = yv | (1 << 63) }
90 ys[n] = yv
91 n = n + 1
92 j = j + 1
93 }
94 // ---- band B: x near 1 (1 +- 2^-k), y ~ 2^(k-12..k-5): dd stress ----
95 j = 0
96 while j < 240 {
97 var r1: i64 = pwg_rng(st)
98 let k: i64 = 20 + (r1 & 15)
99 var xv: i64 = one + ((pwg_rng(st) & 0x000FFFFFFFFFFFFF) >> (k - 12))
100 if (r1 & 64) == 64 { xv = one - ((pwg_rng(st) & 0x000FFFFFFFFFFFFF) >> (k - 12)) }
101 xs[n] = xv
102 let jj: i64 = k - 12 + ((pwg_rng(st) >> 8) & 7)
103 var yv2: i64 = ((1023 + jj) << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
104 if (pwg_rng(st) & 1) == 1 { yv2 = yv2 | (1 << 63) }
105 ys[n] = yv2
106 n = n + 1
107 j = j + 1
108 }
109 // ---- band C: overflow/underflow seams: x in (1, 4), |y| in (512, 2048) ----
110 j = 0
111 while j < 180 {
112 let efx2: i64 = 1023 + (pwg_rng(st) & 1)
113 xs[n] = (efx2 << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
114 let efy2: i64 = 1032 + (pwg_rng(st) & 1)
115 var yv3: i64 = (efy2 << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
116 if (pwg_rng(st) & 1) == 1 { yv3 = yv3 | (1 << 63) }
117 ys[n] = yv3
118 n = n + 1
119 j = j + 1
120 }
121 // ---- band D: wide x (2^-20..2^20), small fractional y ----
122 j = 0
123 while j < 180 {
124 let efx3: i64 = 1003 + (pwg_rng(st) & 31)
125 xs[n] = (efx3 << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
126 let efy3: i64 = 1019 + (pwg_rng(st) & 7)
127 var yv4: i64 = (efy3 << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
128 if (pwg_rng(st) & 1) == 1 { yv4 = yv4 | (1 << 63) }
129 ys[n] = yv4
130 n = n + 1
131 j = j + 1
132 }
133 // ---- band E: integer y (parity paths), x sign random ----
134 j = 0
135 while j < 180 {
136 var xv2: i64 = ((1022 + (pwg_rng(st) & 1)) << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
137 if (pwg_rng(st) & 1) == 1 { xv2 = xv2 | (1 << 63) }
138 xs[n] = xv2
139 var iv: i64 = pwg_rng(st) % 100
140 if iv < 0 { iv = 0 - iv }
141 if (pwg_rng(st) & 1) == 1 { iv = 0 - iv }
142 ys[n] = nx_i64_to_f64(iv)
143 n = n + 1
144 j = j + 1
145 }
146 // ---- band F: subnormal x ----
147 j = 0
148 while j < 60 {
149 var sx: i64 = pwg_rng(st) & 0x000FFFFFFFFFFFFF
150 if sx == 0 { sx = 1 }
151 xs[n] = sx
152 let efy4: i64 = 1020 + (pwg_rng(st) & 3)
153 var yv5: i64 = (efy4 << 52) | (pwg_rng(st) & 0x000FFFFFFFFFFFFF)
154 if (pwg_rng(st) & 1) == 1 { yv5 = yv5 | (1 << 63) }
155 ys[n] = yv5
156 n = n + 1
157 j = j + 1
158 }
159
160 var p_exact: i64 = 0
161 var p_u1: i64 = 0
162 var p_worse: i64 = 0
163 var p_max: i64 = 0
164 var i: i64 = 0
165 while i < n {
166 let pw: i64 = bf_pow_f64(xs[i], ys[i])
167 let pg: i64 = nx_f64_pow(xs[i], ys[i])
168 let pu: i64 = pwg_ulp(pg, pw)
169 if pu == 0 { p_exact = p_exact + 1 }
170 if pu == 1 { p_u1 = p_u1 + 1 }
171 if pu > 1 {
172 p_worse = p_worse + 1
173 if p_worse <= 20 {
174 pwg_puts("PWBAD x=" as *u8); pwg_puthex(xs[i])
175 pwg_puts(" y=" as *u8); pwg_puthex(ys[i])
176 pwg_puts(" exp=" as *u8); pwg_puthex(pw)
177 pwg_puts(" got=" as *u8); pwg_puthex(pg)
178 pwg_puts(" ulp=" as *u8); pwg_putn(pu)
179 pwg_puts("\n" as *u8)
180 }
181 }
182 if pu > p_max { p_max = pu }
183 i = i + 1
184 }
185
186 pwg_puts("POW-ULP total=" as *u8); pwg_putn(n)
187 pwg_puts(" exact=" as *u8); pwg_putn(p_exact)
188 pwg_puts(" ulp1=" as *u8); pwg_putn(p_u1)
189 pwg_puts(" worse=" as *u8); pwg_putn(p_worse)
190 pwg_puts(" max=" as *u8); pwg_putn(p_max)
191 pwg_puts("\n" as *u8)
192 if p_worse == 0 {
193 pwg_puts("POW-GATE verdict=GREEN\n" as *u8)
194 return 0
195 }
196 pwg_puts("POW-GATE verdict=RED\n" as *u8)
197 var rc: i64 = p_worse
198 if rc > 100 { rc = 100 }
199 return rc
200}