_pe_f64pow.nx source
↩ module page · 170 lines · 7090 B
1// AUTHORED BY THE NISHI BUILDER (pattern: MATH_KERNEL / POW) -- no Claude core logic.
2// All constants from the SOVEREIGN bigfloat spec (nx_mathspec_pow). Total domain.
3// dd structure (fdlibm e_pow shape) with exact-rational Taylor families.
4import "nx_syscalls.nx"
5import "nx_tier.nx"
6import "nx_f64.nx"
7import "nx_f64_div.nx"
8import "nx_f64_cvt.nx"
9const PW_HMASK: i64 = 0xFFFFFFFF00000000
10func _pw_ycls(ay: i64) -> i64 {
11 let e: i64 = ((ay >> 52) & 0x7FF) - 1023
12 if e < 0 { return 0 }
13 if e >= 53 { return 2 }
14 let mant: i64 = (ay & 0x000FFFFFFFFFFFFF) | 0x0010000000000000
15 let low: i64 = 52 - e
16 if low > 0 {
17 if (mant & ((1 << low) - 1)) != 0 { return 0 }
18 }
19 if ((mant >> low) & 1) == 1 { return 1 }
20 return 2
21}
22func _pw_flt(a: i64, b: i64) -> i64 {
23 var oa: i64 = a
24 if oa < 0 { oa = (1 << 63) - oa }
25 var ob: i64 = b
26 if ob < 0 { ob = (1 << 63) - ob }
27 if oa < ob { return 1 }
28 return 0
29}
30func _pw_log2m(m_in: i64, n: i64, out: *i64) -> i64 {
31 var m: i64 = m_in
32 var nn: i64 = n
33 var bp: i64 = 4607182418800017408
34 var dph: i64 = 0
35 var dpl: i64 = 0
36 if m < 4608194579719069999 {
37 bp = 4607182418800017408
38 } else {
39 if m < 4610479282544200875 {
40 bp = 4609434218613702656; dph = 4603444092150480896; dpl = 4504109657760791808
41 } else {
42 nn = nn + 1
43 m = m - (1 << 52)
44 }
45 }
46 let u: i64 = nx_f64_sub(m, bp)
47 let v: i64 = nx_f64_div(4607182418800017408, nx_f64_add(m, bp))
48 let ss: i64 = nx_f64_mul(u, v)
49 let sh: i64 = ss & PW_HMASK
50 let th: i64 = nx_f64_add(m, bp) & PW_HMASK
51 let tl: i64 = nx_f64_sub(m, nx_f64_sub(th, bp))
52 let sl: i64 = nx_f64_mul(v, nx_f64_sub(nx_f64_sub(u, nx_f64_mul(sh, th)), nx_f64_mul(sh, tl)))
53 let s2: i64 = nx_f64_mul(ss, ss)
54 var r: i64 = 4593867428597356811
55 r = nx_f64_add(nx_f64_mul(r, s2), 4594314991293244562)
56 r = nx_f64_add(nx_f64_mul(r, s2), 4594856777714582366)
57 r = nx_f64_add(nx_f64_mul(r, s2), 4595526043293882007)
58 r = nx_f64_add(nx_f64_mul(r, s2), 4596373779694328218)
59 r = nx_f64_add(nx_f64_mul(r, s2), 4597482358064142494)
60 r = nx_f64_add(nx_f64_mul(r, s2), 4598584637693219188)
61 r = nx_f64_add(nx_f64_mul(r, s2), 4599676419421066581)
62 r = nx_f64_add(nx_f64_mul(r, s2), 4601392076421969627)
63 r = nx_f64_add(nx_f64_mul(r, s2), 4603579539098121011)
64 r = nx_f64_mul(nx_f64_mul(s2, s2), r)
65 r = nx_f64_add(r, nx_f64_mul(sl, nx_f64_add(sh, ss)))
66 let s2h: i64 = nx_f64_mul(sh, sh)
67 let th3: i64 = nx_f64_add(nx_f64_add(4613937818241073152, s2h), r) & PW_HMASK
68 let tl3: i64 = nx_f64_sub(r, nx_f64_sub(nx_f64_sub(th3, 4613937818241073152), s2h))
69 let u2: i64 = nx_f64_mul(sh, th3)
70 let v2: i64 = nx_f64_add(nx_f64_mul(sl, th3), nx_f64_mul(tl3, ss))
71 let ph: i64 = nx_f64_add(u2, v2) & PW_HMASK
72 let pl: i64 = nx_f64_sub(v2, nx_f64_sub(ph, u2))
73 let zh: i64 = nx_f64_mul(4606838310315229184, ph)
74 let zl: i64 = nx_f64_add(nx_f64_add(nx_f64_mul(4511348162831487992, ph), nx_f64_mul(pl, 4606838314010018813)), dpl)
75 let t: i64 = nx_i64_to_f64(nn)
76 let t1: i64 = nx_f64_add(nx_f64_add(nx_f64_add(zh, zl), dph), t) & PW_HMASK
77 let t2: i64 = nx_f64_sub(zl, nx_f64_sub(nx_f64_sub(nx_f64_sub(t1, t), dph), zh))
78 out[0] = t1
79 out[1] = t2
80 return 0
81}
82func nx_f64_pow(x: i64, y: i64) -> i64 {
83 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
84 let ay: i64 = y & 0x7FFFFFFFFFFFFFFF
85 let xs: i64 = (x >> 63) & 1
86 let ys: i64 = (y >> 63) & 1
87 let one_raw: i64 = 4607182418800017408
88 if ay == 0 { return one_raw }
89 if x == one_raw { return one_raw }
90 var isnan: i64 = 0
91 if ax > 0x7FF0000000000000 { isnan = 1 }
92 if ay > 0x7FF0000000000000 { isnan = 1 }
93 if isnan == 1 { return 0x7FF8000000000000 }
94 if ay == 0x7FF0000000000000 {
95 if ax == one_raw { return one_raw }
96 var big: i64 = 0
97 if ax > one_raw { big = 1 - ys } else { big = ys }
98 if big == 1 { return 0x7FF0000000000000 }
99 return 0
100 }
101 let ycls: i64 = _pw_ycls(ay)
102 if ax == 0x7FF0000000000000 {
103 var rs: i64 = 0
104 if xs == 1 { if ycls == 1 { rs = 1 } }
105 if ys == 0 { return (rs << 63) | 0x7FF0000000000000 }
106 return rs << 63
107 }
108 if ax == 0 {
109 var rs0: i64 = 0
110 if xs == 1 { if ycls == 1 { rs0 = 1 } }
111 if ys == 0 { return rs0 << 63 }
112 return (rs0 << 63) | 0x7FF0000000000000
113 }
114 var rsign: i64 = 0
115 if xs == 1 {
116 if ycls == 0 { return 0x7FF8000000000000 }
117 if ycls == 1 { rsign = 1 }
118 }
119 if y == one_raw { return x }
120 var sig: i64 = ax & 0x000FFFFFFFFFFFFF
121 var ef: i64 = (ax >> 52) & 0x7FF
122 if ef == 0 {
123 ef = 1
124 while sig < 0x0010000000000000 { sig = sig << 1; ef = ef - 1 }
125 sig = sig & 0x000FFFFFFFFFFFFF
126 }
127 let n: i64 = ef - 1023
128 let mraw: i64 = (1023 << 52) | sig
129 let lo2: *i64 = sys_mmap(16) as *i64
130 _pw_log2m(mraw, n, lo2)
131 let t1: i64 = lo2[0]
132 let t2: i64 = lo2[1]
133 let y1: i64 = y & PW_HMASK
134 let pl0: i64 = nx_f64_add(nx_f64_mul(nx_f64_sub(y, y1), t1), nx_f64_mul(y, t2))
135 var ph0: i64 = nx_f64_mul(y1, t1)
136 let z: i64 = nx_f64_add(pl0, ph0)
137 if _pw_flt(z, 4652222813120233472) == 0 { return (rsign << 63) | 0x7FF0000000000000 }
138 if _pw_flt(z, 4652447113492299776 | (1 << 63)) == 1 { return rsign << 63 }
139 var n2: i64 = 0
140 if (z & 0x7FFFFFFFFFFFFFFF) > 4602678819172646912 {
141 var zr: i64 = 0
142 if (z >> 63) == 0 { zr = nx_f64_add(z, 4602678819172646912) } else { zr = nx_f64_sub(z, 4602678819172646912) }
143 n2 = nx_f64_to_i64(zr)
144 ph0 = nx_f64_sub(ph0, nx_i64_to_f64(n2))
145 }
146 let tfr: i64 = nx_f64_add(pl0, ph0) & PW_HMASK
147 let u3: i64 = nx_f64_mul(tfr, 4604418530035630080)
148 let v3: i64 = nx_f64_add(nx_f64_mul(nx_f64_sub(pl0, nx_f64_sub(tfr, ph0)), 4604418534313441775), nx_f64_mul(tfr, 4512570848722726696))
149 let z3: i64 = nx_f64_add(u3, v3)
150 let w3: i64 = nx_f64_sub(v3, nx_f64_sub(z3, u3))
151 let t3: i64 = nx_f64_mul(z3, z3)
152 var pp: i64 = (0 - 4798626848869608100)
153 pp = nx_f64_add(nx_f64_mul(pp, t3), 4448832621486575970)
154 pp = nx_f64_add(nx_f64_mul(pp, t3), (0 - 4750690651387713921))
155 pp = nx_f64_add(nx_f64_mul(pp, t3), 4496398441126497982)
156 pp = nx_f64_add(nx_f64_mul(pp, t3), (0 - 4702957064589873397))
157 pp = nx_f64_add(nx_f64_mul(pp, t3), 4544508515414250855)
158 pp = nx_f64_add(nx_f64_mul(pp, t3), (0 - 4654820494858425321))
159 pp = nx_f64_add(nx_f64_mul(pp, t3), 4595172819793696085)
160 let t1p: i64 = nx_f64_sub(z3, nx_f64_mul(t3, pp))
161 let r3: i64 = nx_f64_sub(nx_f64_div(nx_f64_mul(z3, t1p), nx_f64_sub(t1p, 4611686018427387904)), w3)
162 let zf: i64 = nx_f64_sub(4607182418800017408, nx_f64_sub(r3, z3))
163 let efz: i64 = (zf >> 52) & 0x7FF
164 let eo: i64 = efz + n2
165 if eo >= 2047 { return (rsign << 63) | 0x7FF0000000000000 }
166 if eo >= 1 { return (rsign << 63) | (zf + (n2 << 52)) }
167 let m1: i64 = zf - (900 << 52)
168 let pw2: i64 = (n2 + 1923) << 52
169 return (rsign << 63) | nx_f64_mul(m1, pw2)
170}