code wiki / _hdl_build / nx_mathspec_pow.nx
nx_mathspec_pow.nx source
↩ module page · 214 lines · 7176 B
1// nx_mathspec_pow.nx -- SOVEREIGN SPEC EMITTER for the POW kernel: every constant
2// computed on the 120-bit bigfloat; written to _pm_pow_spec.nx for the
3// MATH_KERNEL/POW emitter. No python, no minimax tables: the kernel's
4// correction polynomials use EXACT-RATIONAL Taylor families
5// L_j = 3/(2j+3) (log2 atanh-correction, |s| <= 0.101 -> 10 terms)
6// P_k = 2*B_{2k}/(2k)! (exp R(z)=z(e^z+1)/(e^z-1) Bernoulli form, 8 terms)
7// where B_2..B_16 = 1/6 -1/30 1/42 -1/30 5/66 -691/2730 7/6 -3617/510 (DLMF 24.2,
8// mathematical definitions like 1/k!). (2k)! division runs as a loop of small
9// divisors 2..2k (bf_div_small bound). The sqrt(3/2)/sqrt(3) anchor cuts come
10// from 53-step bigfloat binary search on sig^2 vs 3*2^103 / 3*2^104.
11// 32-bit hi/lo splits (dp, cp, lg2) mask AFTER rounding; remainder via bf_sub.
12//
13// Layout: c[0]=one c[1]=half c[2]=two c[3]=cut_sqrt15 c[4]=cut_sqrt3 c[5]=bp1
14// c[6]=dp_h c[7]=dp_l (log2(1.5) split) c[8]=cp_h c[9]=cp_l c[10]=cp(2/(3ln2))
15// c[11]=lg2 c[12]=lg2_h c[13]=lg2_l c[14..23]=L10..L1 (hi-first)
16// c[24..31]=P8..P1 (hi-first, signed) c[32]=three c[33]=1025.0 c[34]=1076.0
17// license_tier: ORIGINAL
18import "nx_syscalls.nx"
19import "nx_pattern_emit.nx"
20import "nx_bigfloat120.nx"
21import "nx_bigfloat120_div.nx"
22import "nx_bigfloat120_exp.nx"
23const K_MAGIC_2730: i64 = 2730
24const K_MAGIC_3617: i64 = 3617
25const K_MAGIC_1025: i64 = 1025
26const K_MAGIC_1076: i64 = 1076
27
28func mpw_wlit(fd: i64, v: i64) -> i64 {
29 if v < 0 {
30 pe_w(fd, "0 - " as *u8)
31 pe_wn(fd, 0 - v)
32 return 0
33 }
34 pe_wn(fd, v)
35 return 0
36}
37
38func mpw_emit_row(fd: i64, idx: i64, v: i64) -> i64 {
39 pe_w(fd, " c[" as *u8); pe_wn(fd, idx); pe_w(fd, "] = " as *u8)
40 mpw_wlit(fd, v)
41 pe_w(fd, "\n" as *u8)
42 return 0
43}
44
45// 32-bit split rows: hi = fl(v) with low 32 raw bits cleared (toward zero, so
46// the remainder stays positive on the positive-only core), lo = fl(v - hi).
47func mpw_emit_split(fd: i64, idxhi: i64, idxlo: i64, v: *i64) -> i64 {
48 let hi: i64 = bf_to_f64(v, 0, 0) & 0xFFFFFFFF00000000
49 let hib: *i64 = bf_new()
50 bf_set_f64(hib, hi)
51 let rem: *i64 = bf_new()
52 bf_sub(rem, v, hib) // v >= hib by truncation
53 mpw_emit_row(fd, idxhi, hi)
54 mpw_emit_row(fd, idxlo, bf_to_f64(rem, 0, 0))
55 return 0
56}
57
58// smallest 53-bit sig with sig*sig >= rad (rad = small*2^eshift as bigfloat);
59// returned as an f64 raw in the [1,2) binade: (1023<<52) | (sig - 2^52).
60func mpw_sqrt_cut(small: i64, eshift: i64) -> i64 {
61 let rad: *i64 = bf_new()
62 bf_set_int(rad, small)
63 rad[0] = rad[0] + eshift
64 let s: *i64 = bf_new()
65 let sq: *i64 = bf_new()
66 var lo: i64 = 1 << 52
67 var hi: i64 = 1 << 53
68 while lo < hi {
69 let mid: i64 = lo + ((hi - lo) >> 1)
70 bf_set_int(s, mid)
71 bf_mul(sq, s, s)
72 if bf_cmp(sq, rad) >= 0 { hi = mid } else { lo = mid + 1 }
73 }
74 return (1023 << 52) | (lo - (1 << 52))
75}
76
77// atanh(1/5) * 2 = ln(1.5), then / ln2 -> log2(1.5)
78func mpw_log2_15(out: *i64, ln2: *i64) -> i64 {
79 let one: *i64 = bf_new()
80 bf_set_int(one, 1)
81 let x: *i64 = bf_new()
82 bf_div_small(x, one, 5)
83 let z: *i64 = bf_new()
84 bf_mul(z, x, x)
85 let term: *i64 = bf_new()
86 bf_copy(term, x)
87 let sum: *i64 = bf_new()
88 bf_copy(sum, x)
89 let tmul: *i64 = bf_new()
90 let tdiv: *i64 = bf_new()
91 let snew: *i64 = bf_new()
92 var k: i64 = 1
93 while k <= 30 {
94 bf_mul(tmul, term, z)
95 bf_copy(term, tmul)
96 bf_div_small(tdiv, term, 2 * k + 1)
97 bf_add(snew, sum, tdiv)
98 bf_copy(sum, snew)
99 k = k + 1
100 }
101 sum[0] = sum[0] + 1 // * 2
102 return bf_div(out, sum, ln2)
103}
104
105// P_k = 2*|B_2k| / (2k)! with sign (-1)^(k+1); divides by 2..2k one factor at
106// a time (each < bf_div_small's bound).
107func mpw_pk(out: *i64, bnum: i64, bden: i64, twok: i64) -> i64 {
108 let t: *i64 = bf_new()
109 bf_set_int(t, 2 * bnum)
110 let u: *i64 = bf_new()
111 bf_div_small(u, t, bden)
112 var d: i64 = 2
113 while d <= twok {
114 bf_div_small(t, u, d)
115 bf_copy(u, t)
116 d = d + 1
117 }
118 return bf_copy(out, u)
119}
120
121func main() -> i64 {
122 let fd: i64 = sys_openat_wr("runtime/_hdl_build/_pm_pow_spec.nx" as *u8, 0x1a4)
123 if fd < 0 { sys_exit(2) }
124 pe_w(fd, "// _pm_pow_spec.nx -- GENERATED by nx_mathspec_pow (SOVEREIGN bigfloat,\n" as *u8)
125 pe_w(fd, "// no python). DO NOT HAND-EDIT. Layout: see nx_mathspec_pow.nx header.\n" as *u8)
126 pe_w(fd, "import \"nx_syscalls.nx\"\n\nfunc pm_pow_spec_fill(c: *i64) -> i64 {\n" as *u8)
127
128 let one: *i64 = bf_new()
129 bf_set_int(one, 1)
130 mpw_emit_row(fd, 0, bf_to_f64(one, 0, 0))
131 let half: *i64 = bf_new()
132 bf_copy(half, one)
133 half[0] = half[0] - 1
134 mpw_emit_row(fd, 1, bf_to_f64(half, 0, 0))
135 let two: *i64 = bf_new()
136 bf_set_int(two, 2)
137 mpw_emit_row(fd, 2, bf_to_f64(two, 0, 0))
138
139 // anchor cuts: sqrt(3/2) = sqrt(3*2^103)/2^52, sqrt(3) = sqrt(3*2^104)/2^52
140 mpw_emit_row(fd, 3, mpw_sqrt_cut(3, 103))
141 mpw_emit_row(fd, 4, mpw_sqrt_cut(3, 104))
142
143 let bp1: *i64 = bf_new()
144 bf_set_int(bp1, 3)
145 bp1[0] = bp1[0] - 1 // 1.5
146 mpw_emit_row(fd, 5, bf_to_f64(bp1, 0, 0))
147
148 let ln2: *i64 = bf_new()
149 bf_ln2(ln2)
150 let dp: *i64 = bf_new()
151 mpw_log2_15(dp, ln2)
152 mpw_emit_split(fd, 6, 7, dp)
153
154 // cp = 2/(3*ln2)
155 let t3: *i64 = bf_new()
156 bf_div_small(t3, two, 3)
157 let cp: *i64 = bf_new()
158 bf_div(cp, t3, ln2)
159 mpw_emit_split(fd, 8, 9, cp)
160 mpw_emit_row(fd, 10, bf_to_f64(cp, 0, 0))
161
162 mpw_emit_row(fd, 11, bf_to_f64(ln2, 0, 0))
163 mpw_emit_split(fd, 12, 13, ln2)
164
165 // L10..L1 hi-first: L_j = 3/(2j+3)
166 let three: *i64 = bf_new()
167 bf_set_int(three, 3)
168 var j: i64 = 10
169 var idx: i64 = 14
170 while j >= 1 {
171 let lj: *i64 = bf_new()
172 bf_div_small(lj, three, 2 * j + 3)
173 mpw_emit_row(fd, idx, bf_to_f64(lj, 0, 0))
174 idx = idx + 1
175 j = j - 1
176 }
177
178 // P8..P1 hi-first: |B_2k| num/den pairs, sign (-1)^(k+1)
179 let bn: *i64 = sys_mmap(8 * 10) as *i64
180 let bd: *i64 = sys_mmap(8 * 10) as *i64
181 bn[1] = 1; bd[1] = 6
182 bn[2] = 1; bd[2] = 30
183 bn[3] = 1; bd[3] = 42
184 bn[4] = 1; bd[4] = 30
185 bn[5] = 5; bd[5] = 66
186 bn[6] = 691; bd[6] = K_MAGIC_2730
187 bn[7] = 7; bd[7] = 6
188 bn[8] = K_MAGIC_3617; bd[8] = 510
189 var k: i64 = 8
190 idx = 24
191 while k >= 1 {
192 let pk: *i64 = bf_new()
193 mpw_pk(pk, bn[k], bd[k], 2 * k)
194 var bits: i64 = bf_to_f64(pk, 0, 0)
195 if (k & 1) == 0 { bits = bits | (1 << 63) } // sign (-1)^(k+1)
196 mpw_emit_row(fd, idx, bits)
197 idx = idx + 1
198 k = k - 1
199 }
200
201 mpw_emit_row(fd, 32, bf_to_f64(three, 0, 0))
202 let b25: *i64 = bf_new()
203 bf_set_int(b25, K_MAGIC_1025)
204 mpw_emit_row(fd, 33, bf_to_f64(b25, 0, 0))
205 let b76: *i64 = bf_new()
206 bf_set_int(b76, K_MAGIC_1076)
207 mpw_emit_row(fd, 34, bf_to_f64(b76, 0, 0))
208
209 pe_w(fd, " return 35\n}\n" as *u8)
210 sys_close(fd)
211 pe_w(1, "SPEC EMITTED: _pm_pow_spec.nx (35 consts, all from sovereign bigfloat)\n" as *u8)
212 sys_exit(0)
213 return 0
214}