code wiki / _hdl_build / nx_mathspec_sinhcosh.nx
nx_mathspec_sinhcosh.nx source
↩ module page · 156 lines · 4749 B
1// nx_mathspec_sinhcosh.nx -- SOVEREIGN SPEC EMITTER for the SINH/COSH kernels:
2// every constant from the 120-bit bigfloat. No python, no minimax:
3// exp-dd core reuses the POW families (ln2 32-bit hi/lo split, Bernoulli
4// P_k = 2*B_2k/(2k)!) plus INV_LN2 for the reduction multiply;
5// Taylor paths use exact factorial reciprocals 1/(2k+1)! (sinh) and
6// 1/(2k)! (cosh), divided factor-by-factor under the bf_div_small bound.
7//
8// Layout: c[0]=one c[1]=half c[2]=two c[3]=INV_LN2 c[4]=ln2_h c[5]=ln2_l
9// c[6..13]=P8..P1 (hi-first, signed)
10// c[14..23]=sinh coeffs hi-first: 1/(2k+1)! k=10..1
11// c[24..33]=cosh coeffs hi-first: 1/(2k)! k=10..1
12// c[34]=28.0 (e^-2x sub-ulp cut) c[35]=711.0 (hard overflow cut)
13// license_tier: ORIGINAL
14import "nx_syscalls.nx"
15import "nx_pattern_emit.nx"
16import "nx_bigfloat120.nx"
17import "nx_bigfloat120_div.nx"
18import "nx_bigfloat120_exp.nx"
19const K_MAGIC_2730: i64 = 2730
20const K_MAGIC_3617: i64 = 3617
21
22func msh_wlit(fd: i64, v: i64) -> i64 {
23 if v < 0 {
24 pe_w(fd, "0 - " as *u8)
25 pe_wn(fd, 0 - v)
26 return 0
27 }
28 pe_wn(fd, v)
29 return 0
30}
31
32func msh_emit_row(fd: i64, idx: i64, v: i64) -> i64 {
33 pe_w(fd, " c[" as *u8); pe_wn(fd, idx); pe_w(fd, "] = " as *u8)
34 msh_wlit(fd, v)
35 pe_w(fd, "\n" as *u8)
36 return 0
37}
38
39// P_k = 2*|B_2k|/(2k)! with sign (-1)^(k+1)
40func msh_pk(out: *i64, bnum: i64, bden: i64, twok: i64) -> i64 {
41 let t: *i64 = bf_new()
42 bf_set_int(t, 2 * bnum)
43 let u: *i64 = bf_new()
44 bf_div_small(u, t, bden)
45 var d: i64 = 2
46 while d <= twok {
47 bf_div_small(t, u, d)
48 bf_copy(u, t)
49 d = d + 1
50 }
51 return bf_copy(out, u)
52}
53
54// 1/m! via factor-by-factor division
55func msh_invfact(out: *i64, m: i64) -> i64 {
56 let t: *i64 = bf_new()
57 bf_set_int(t, 1)
58 let u: *i64 = bf_new()
59 bf_copy(u, t)
60 var d: i64 = 2
61 while d <= m {
62 bf_div_small(t, u, d)
63 bf_copy(u, t)
64 d = d + 1
65 }
66 return bf_copy(out, u)
67}
68
69func main() -> i64 {
70 let fd: i64 = sys_openat_wr("runtime/_hdl_build/_pm_sinhcosh_spec.nx" as *u8, 0x1a4)
71 if fd < 0 { sys_exit(2) }
72 pe_w(fd, "// _pm_sinhcosh_spec.nx -- GENERATED by nx_mathspec_sinhcosh (SOVEREIGN\n" as *u8)
73 pe_w(fd, "// bigfloat, no python). DO NOT HAND-EDIT. Layout: see emitter header.\n" as *u8)
74 pe_w(fd, "import \"nx_syscalls.nx\"\n\nfunc pm_sinhcosh_spec_fill(c: *i64) -> i64 {\n" as *u8)
75
76 let one: *i64 = bf_new()
77 bf_set_int(one, 1)
78 msh_emit_row(fd, 0, bf_to_f64(one, 0, 0))
79 let half: *i64 = bf_new()
80 bf_copy(half, one)
81 half[0] = half[0] - 1
82 msh_emit_row(fd, 1, bf_to_f64(half, 0, 0))
83 let two: *i64 = bf_new()
84 bf_set_int(two, 2)
85 msh_emit_row(fd, 2, bf_to_f64(two, 0, 0))
86
87 let ln2: *i64 = bf_new()
88 bf_ln2(ln2)
89 let inv: *i64 = bf_new()
90 bf_div(inv, one, ln2)
91 msh_emit_row(fd, 3, bf_to_f64(inv, 0, 0))
92 let lh: i64 = bf_to_f64(ln2, 0, 0) & 0xFFFFFFFF00000000
93 let lhb: *i64 = bf_new()
94 bf_set_f64(lhb, lh)
95 let lrem: *i64 = bf_new()
96 bf_sub(lrem, ln2, lhb)
97 msh_emit_row(fd, 4, lh)
98 msh_emit_row(fd, 5, bf_to_f64(lrem, 0, 0))
99
100 // P8..P1 hi-first (Bernoulli |B_2k| num/den, sign (-1)^(k+1))
101 let bn: *i64 = sys_mmap(8 * 10) as *i64
102 let bd: *i64 = sys_mmap(8 * 10) as *i64
103 bn[1] = 1; bd[1] = 6
104 bn[2] = 1; bd[2] = 30
105 bn[3] = 1; bd[3] = 42
106 bn[4] = 1; bd[4] = 30
107 bn[5] = 5; bd[5] = 66
108 bn[6] = 691; bd[6] = K_MAGIC_2730
109 bn[7] = 7; bd[7] = 6
110 bn[8] = K_MAGIC_3617; bd[8] = 510
111 var k: i64 = 8
112 var idx: i64 = 6
113 while k >= 1 {
114 let pk: *i64 = bf_new()
115 msh_pk(pk, bn[k], bd[k], 2 * k)
116 var bits: i64 = bf_to_f64(pk, 0, 0)
117 if (k & 1) == 0 { bits = bits | (1 << 63) }
118 msh_emit_row(fd, idx, bits)
119 idx = idx + 1
120 k = k - 1
121 }
122
123 // sinh coeffs hi-first: 1/(2k+1)! k=10..1
124 k = 10
125 idx = 14
126 while k >= 1 {
127 let sc: *i64 = bf_new()
128 msh_invfact(sc, 2 * k + 1)
129 msh_emit_row(fd, idx, bf_to_f64(sc, 0, 0))
130 idx = idx + 1
131 k = k - 1
132 }
133 // cosh coeffs hi-first: 1/(2k)! k=10..1
134 k = 10
135 idx = 24
136 while k >= 1 {
137 let cc: *i64 = bf_new()
138 msh_invfact(cc, 2 * k)
139 msh_emit_row(fd, idx, bf_to_f64(cc, 0, 0))
140 idx = idx + 1
141 k = k - 1
142 }
143
144 let c28: *i64 = bf_new()
145 bf_set_int(c28, 28)
146 msh_emit_row(fd, 34, bf_to_f64(c28, 0, 0))
147 let c711: *i64 = bf_new()
148 bf_set_int(c711, 711)
149 msh_emit_row(fd, 35, bf_to_f64(c711, 0, 0))
150
151 pe_w(fd, " return 36\n}\n" as *u8)
152 sys_close(fd)
153 pe_w(1, "SPEC EMITTED: _pm_sinhcosh_spec.nx (36 consts, all from sovereign bigfloat)\n" as *u8)
154 sys_exit(0)
155 return 0
156}