code wiki / _hdl_build / nx_mathspec_atan.nx
`nx_mathspec_atan.nx` source

↩ module page · 160 lines · 5153 B
1// nx_mathspec_atan.nx -- SOVEREIGN SPEC EMITTER for the ATAN kernel: every
2// constant computed on the 120-bit bigfloat (oracle organ nx_bigfloat120_atan),
3// written to _pm_atan_spec.nx for the MATH_KERNEL/ATAN emitter.  No python.
4//
5// Kernel scheme (accurate-table): anchors c_j = j/8, j = 0..20.  For
6// x < 41/16 the kernel takes j = nearest(8x), t = (x - c_j)/(1 + x*c_j)
7// (|t| <= 1/16, x - c_j Sterbenz-exact), atan(x) = A_hi[j] + (A_lo[j] + P(t))
8// with P an 8-coeff odd Taylor poly (truncation < 2^-64 rel).  For x >= 41/16
9// it folds atan(x) = pi/2 - atan(1/x); there 1/x < 0.4 so the fold target sits
10// 2 binades below the result and the transferred error is damped 4x.
11//
12// Layout: c[0]=one c[1]=half c[2]=eight c[3]=pi2_hi c[4]=pi2_lo c[5]=cut(41/16)
13//         c[6..26]   anchors j/8 (j=0..20)
14//         c[27..47]  A_hi[j] = fl(atan(j/8))
15//         c[48..68]  A_lo[j] = fl(atan(j/8) - A_hi[j])   (signed)
16//         c[69..76]  poly coeffs hi-first: -1/15 1/13 -1/11 1/9 -1/7 1/5 -1/3 1
17// NEGATIVE literals emitted in `0 - N` form (unary-minus literal landmine).
18// license_tier: ORIGINAL
19import "nx_syscalls.nx"
20import "nx_pattern_emit.nx"
21import "nx_bigfloat120.nx"
22import "nx_bigfloat120_div.nx"
23import "nx_bigfloat120_trig.nx"
24import "nx_bigfloat120_atan.nx"
25const K_MAGIC_999999: i64 = 999999
26const K_MAGIC_5625: i64 = 5625
27
28func mat_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 mat_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    mat_wlit(fd, v)
41    pe_w(fd, "\n" as *u8)
42    return 0
43}
44
45// atan(j/8) on the bigfloat for j = 0..20 (j > 8 folds through pi/2 - atan(8/j))
46func mat_atan_j8(out: *i64, j: i64) -> i64 {
47    if j == 0 { out[0] = 0 - K_MAGIC_999999; out[1] = 0; out[2] = 0; return 0 }
48    let one: *i64 = bf_new()
49    bf_set_int(one, 1)
50    if j <= 8 {
51        let x: *i64 = bf_new()
52        bf_set_int(x, j)
53        x[0] = x[0] - 3
54        return bf_atan01(out, x)
55    }
56    let eightb: *i64 = bf_new()
57    bf_set_int(eightb, 8)
58    let inv: *i64 = bf_new()
59    bf_div_small(inv, eightb, j)                 // 8/j < 1
60    let s: *i64 = bf_new()
61    bf_atan01(s, inv)
62    let p2: *i64 = bf_new()
63    bf_pi(p2)
64    p2[0] = p2[0] - 1
65    return bf_sub(out, p2, s)
66}
67
68// split a positive bigfloat into f64 hi + signed f64 lo; writes both rows
69func mat_emit_hilo(fd: i64, idxhi: i64, idxlo: i64, a: *i64) -> i64 {
70    if bf_is_zero(a) == 1 {
71        mat_emit_row(fd, idxhi, 0)
72        mat_emit_row(fd, idxlo, 0)
73        return 0
74    }
75    let hi: i64 = bf_to_f64(a, 0, 0)
76    let hib: *i64 = bf_new()
77    bf_set_f64(hib, hi)
78    let rem: *i64 = bf_new()
79    var losig: i64 = 0
80    if bf_cmp(a, hib) >= 0 {
81        bf_sub(rem, a, hib)
82    } else {
83        bf_sub(rem, hib, a)
84        losig = 1
85    }
86    mat_emit_row(fd, idxhi, hi)
87    mat_emit_row(fd, idxlo, bf_to_f64(rem, losig, 0))
88    return 0
89}
90
91func main() -> i64 {
92    let fd: i64 = sys_openat_wr("runtime/_hdl_build/_pm_atan_spec.nx" as *u8, 0x1a4)
93    if fd < 0 { sys_exit(2) }
94    pe_w(fd, "// _pm_atan_spec.nx -- GENERATED by nx_mathspec_atan (SOVEREIGN bigfloat,\n" as *u8)
95    pe_w(fd, "// no python). DO NOT HAND-EDIT. Layout: see nx_mathspec_atan.nx header.\n" as *u8)
96    pe_w(fd, "import \"nx_syscalls.nx\"\n\nfunc pm_atan_spec_fill(c: *i64) -> i64 {\n" as *u8)
97
98    let one: *i64 = bf_new()
99    bf_set_int(one, 1)
100    mat_emit_row(fd, 0, bf_to_f64(one, 0, 0))
101    let half: *i64 = bf_new()
102    bf_copy(half, one)
103    half[0] = half[0] - 1
104    mat_emit_row(fd, 1, bf_to_f64(half, 0, 0))
105    let eightb: *i64 = bf_new()
106    bf_set_int(eightb, 8)
107    mat_emit_row(fd, 2, bf_to_f64(eightb, 0, 0))
108
109    let p2: *i64 = bf_new()
110    bf_pi(p2)
111    p2[0] = p2[0] - 1
112    mat_emit_hilo(fd, 3, 4, p2)
113
114    let cut: *i64 = bf_new()
115    bf_set_int(cut, 41)
116    cut[0] = cut[0] - 4                          // 41/16 = 2.K_MAGIC_5625
117    mat_emit_row(fd, 5, bf_to_f64(cut, 0, 0))
118
119    // anchors c_j = j/8 (exact in f64)
120    var j: i64 = 0
121    while j <= 20 {
122        if j == 0 {
123            mat_emit_row(fd, 6, 0)
124        } else {
125            let cj: *i64 = bf_new()
126            bf_set_int(cj, j)
127            cj[0] = cj[0] - 3
128            mat_emit_row(fd, 6 + j, bf_to_f64(cj, 0, 0))
129        }
130        j = j + 1
131    }
132
133    // A_hi / A_lo tables
134    j = 0
135    while j <= 20 {
136        let aj: *i64 = bf_new()
137        mat_atan_j8(aj, j)
138        mat_emit_hilo(fd, 27 + j, 48 + j, aj)
139        j = j + 1
140    }
141
142    // poly coeffs hi-first: k = 7..0, coeff = (-1)^k / (2k+1)
143    var k: i64 = 7
144    var idx: i64 = 69
145    while k >= 0 {
146        let cf: *i64 = bf_new()
147        bf_div_small(cf, one, 2 * k + 1)
148        var bits: i64 = bf_to_f64(cf, 0, 0)
149        if (k & 1) == 1 { bits = bits | (1 << 63) }
150        mat_emit_row(fd, idx, bits)
151        idx = idx + 1
152        k = k - 1
153    }
154
155    pe_w(fd, "    return 77\n}\n" as *u8)
156    sys_close(fd)
157    pe_w(1, "SPEC EMITTED: _pm_atan_spec.nx (77 consts, all from sovereign bigfloat atan)\n" as *u8)
158    sys_exit(0)
159    return 0
160}
nx_mathspec_atan.nx source

`nx_mathspec_atan.nx` source
