code wiki / _hdl_build / nx_mathspec_erf.nx
nx_mathspec_erf.nx source
↩ module page · 150 lines · 5723 B
1// nx_mathspec_erf.nx -- SOVEREIGN SPEC EMITTER for the ERF kernel (DLMF ch7 rung 1,
2// ME2-ERF-001). Every constant from the 120-bit bigfloat -- no python, no minimax:
3// 2/sqrt(pi) via bf_pi (Machin) + NEWTON square root (8 iterations from the f64
4// seed, doubling precision each step: bf_div exists, so s <- (s + pi/s)/2);
5// Taylor reciprocals r_n = 1/(n!(2n+1)) by factor-by-factor bf_div_small with
6// the sign (-1)^n BAKED into the emitted bits (the sinhcosh P_k sign gene).
7// REFUSAL RAILS (refuse > fabricate):
8// R1: (2/sqrt(pi)) * sqrt(pi) must equal 2 to ~100 bits, else exit 2 no emit;
9// R2: the derived 2/sqrt(pi) f64 bits must equal the KNOWN-GOOD libm constant
10// 0x3FF20DD750429B6D bit-exactly (derive sovereignly, verify vs the world).
11// Layout: c[0]=one c[1]=two_over_sqrt_pi c[2]=tosp_hi(32-bit split) c[3]=tosp_lo
12// c[4..35]=32 SIGNED Taylor reciprocals hi-first: (-1)^n/(n!(2n+1)) n=31..0
13// (24 terms left the tail at ~2^-42 at the cut; 32 puts it past 2^-67 -- ULP-clean)
14// c[36]=series_cut 1.75 c[37]=one_cut 6.0 (|x|>=6 -> erf=+-1 sub-ulp)
15// Emits runtime/_hdl_build/_pm_erf_spec.nx :: pm_erf_spec_fill(c)->30
16// license_tier: ORIGINAL
17import "nx_syscalls.nx"
18import "nx_pattern_emit.nx"
19import "nx_bigfloat120.nx"
20import "nx_bigfloat120_div.nx"
21import "nx_bigfloat120_trig.nx"
22
23func mse_wlit(fd: i64, v: i64) -> i64 {
24 if v < 0 {
25 pe_w(fd, "0 - " as *u8)
26 pe_wn(fd, 0 - v)
27 return 0
28 }
29 pe_wn(fd, v)
30 return 0
31}
32
33func mse_emit_row(fd: i64, idx: i64, v: i64) -> i64 {
34 pe_w(fd, " c[" as *u8); pe_wn(fd, idx); pe_w(fd, "] = " as *u8)
35 mse_wlit(fd, v)
36 pe_w(fd, "\n" as *u8)
37 return 0
38}
39
40// r_n = 1/(n! * (2n+1)) by factor-by-factor division (the invfact gene + one more)
41func mse_recip(out: *i64, n: i64) -> i64 {
42 let t: *i64 = bf_new()
43 bf_set_int(t, 1)
44 let u: *i64 = bf_new()
45 bf_copy(u, t)
46 var d: i64 = 2
47 while d <= n {
48 bf_div_small(t, u, d)
49 bf_copy(u, t)
50 d = d + 1
51 }
52 bf_div_small(t, u, 2 * n + 1)
53 return bf_copy(out, t)
54}
55
56// Newton square root: s <- (s + a/s)/2, 8 iterations from an f64 seed
57func mse_sqrt(out: *i64, a: *i64, seed_raw: i64) -> i64 {
58 let s: *i64 = bf_new()
59 bf_set_f64(s, seed_raw)
60 let q: *i64 = bf_new()
61 let sum: *i64 = bf_new()
62 var i: i64 = 0
63 while i < 8 {
64 bf_div(q, a, s)
65 bf_add(sum, s, q)
66 bf_div_small(s, sum, 2)
67 i = i + 1
68 }
69 return bf_copy(out, s)
70}
71
72func main() -> i64 {
73 // ---- derive 2/sqrt(pi) sovereignly ----
74 let pi: *i64 = bf_new()
75 bf_pi(pi)
76 let spi: *i64 = bf_new()
77 // f64 seed for sqrt(pi) ~ 1.7724538509055160 (0x3FFC5BF891B4EF6A); Newton
78 // needs only ~50 good bits -- 8 doublings overshoot 120 comfortably
79 mse_sqrt(spi, pi, 0x3FFC5BF891B4EF6A)
80 let two: *i64 = bf_new()
81 bf_set_int(two, 2)
82 let tosp: *i64 = bf_new()
83 bf_div(tosp, two, spi)
84 // ---- R1: tosp * spi == 2 to ~100 bits ----
85 let chk: *i64 = bf_new()
86 bf_mul(chk, tosp, spi)
87 let dif: *i64 = bf_new()
88 // bigfloats are magnitude-only (signs live in emitted bits) -- order the
89 // subtraction by comparison or the smaller-minus-larger case underflows
90 if bf_cmp(chk, two) >= 0 { bf_sub(dif, chk, two) } else { bf_sub(dif, two, chk) }
91 var db: i64 = bf_to_f64(dif, 0, 0)
92 db = db & 0x7FFFFFFFFFFFFFFF
93 if db >= (923 << 52) {
94 pe_w(1, "REFUSED: Newton sqrt self-check failed (|2/sqrt(pi)*sqrt(pi)-2| too large)\n" as *u8)
95 sys_exit(2)
96 return 2
97 }
98 // ---- R2: derived bits must equal the known-good libm constant ----
99 let tosp_bits: i64 = bf_to_f64(tosp, 0, 0)
100 if tosp_bits != 0x3FF20DD750429B6D {
101 pe_w(1, "REFUSED: derived 2/sqrt(pi) bits differ from known-good 0x3FF20DD750429B6D\n" as *u8)
102 sys_exit(2)
103 return 2
104 }
105 // ---- emit the spec ----
106 let fd: i64 = sys_openat_wr("runtime/_hdl_build/_pm_erf_spec.nx" as *u8, 0x1a4)
107 if fd < 0 { sys_exit(2) }
108 pe_w(fd, "// _pm_erf_spec.nx -- GENERATED by nx_mathspec_erf (SOVEREIGN bigfloat:\n" as *u8)
109 pe_w(fd, "// bf_pi Machin + Newton sqrt + factor-division). DO NOT HAND-EDIT.\n" as *u8)
110 pe_w(fd, "// Layout: see emitter header. Both refusal rails passed at emit time.\n" as *u8)
111 pe_w(fd, "import \"nx_syscalls.nx\"\n\nfunc pm_erf_spec_fill(c: *i64) -> i64 {\n" as *u8)
112 let one: *i64 = bf_new()
113 bf_set_int(one, 1)
114 mse_emit_row(fd, 0, bf_to_f64(one, 0, 0))
115 mse_emit_row(fd, 1, tosp_bits)
116 // dd split of 2/sqrt(pi): hi = top-32 mantissa bits, lo = exact residual
117 let th: i64 = tosp_bits & 0xFFFFFFFF00000000
118 let thb: *i64 = bf_new()
119 bf_set_f64(thb, th)
120 let trem: *i64 = bf_new()
121 bf_sub(trem, tosp, thb)
122 mse_emit_row(fd, 2, th)
123 mse_emit_row(fd, 3, bf_to_f64(trem, 0, 0))
124 // 32 signed Taylor reciprocals hi-first: (-1)^n / (n!(2n+1)), n = 31..0
125 var n: i64 = 31
126 var idx: i64 = 4
127 while n >= 0 {
128 let r: *i64 = bf_new()
129 mse_recip(r, n)
130 var bits: i64 = bf_to_f64(r, 0, 0)
131 if (n & 1) == 1 { bits = bits | (1 << 63) }
132 mse_emit_row(fd, idx, bits)
133 idx = idx + 1
134 n = n - 1
135 }
136 // cuts: series domain 1.75 = 7/4; one-domain 6.0
137 let c7: *i64 = bf_new()
138 bf_set_int(c7, 7)
139 let cq: *i64 = bf_new()
140 bf_div_small(cq, c7, 4)
141 mse_emit_row(fd, 36, bf_to_f64(cq, 0, 0))
142 let c6: *i64 = bf_new()
143 bf_set_int(c6, 6)
144 mse_emit_row(fd, 37, bf_to_f64(c6, 0, 0))
145 pe_w(fd, " return 38\n}\n" as *u8)
146 sys_close(fd)
147 pe_w(1, "SPEC EMITTED: _pm_erf_spec.nx (38 consts; rails R1+R2 held: sovereign == known-good)\n" as *u8)
148 sys_exit(0)
149 return 0
150}