nx_bigfloat120_pow.nx source
↩ module page · 239 lines · 8047 B
1// nx_bigfloat120_pow.nx -- SOVEREIGN pow oracle: binary64 pair in, binary64 out,
2// pow(x,y) = exp(y * ln|x|) entirely on the 120-bit bigfloat (sign by the
3// odd-integer-y rule). Margin: |y*ln x| <= 1100 for any finite nonzero result,
4// so the bigfloat's 2^-110 relative error gives the exp argument an ABSOLUTE
5// error < 2^-99 -- the result is correct to ~2^-99 relative, 45 bits past the
6// f64 half-ulp decision boundary. |v| > 1100 is decided as hard inf/zero
7// (the finite-result band ends at |ln r| ~ 745.2); inside the band bf_to_f64
8// handles overflow -> inf and gradual underflow -> subnormal/zero exactly.
9//
10// The IEEE-754/C99 special matrix is decided in pure bit logic BEFORE any
11// bigfloat runs (pow(x,+-0)=1 and pow(+1,y)=1 even for NaN partners).
12//
13// Self-anchor (see _bf_pow_gate_authored): pow(x,2)=fl(x*x), pow(x,-1)=fl(1/x),
14// pow(x,1/2)=sqrt(x) against the BLESSED f64 mul/div/sqrt organs (correctly-
15// rounded equalities), exact power-of-two ladders incl. the subnormal floor and
16// the overflow ceiling, plus the full special matrix.
17// license_tier: ORIGINAL
18
19import "nx_syscalls.nx"
20import "nx_tier.nx"
21import "nx_bigfloat120.nx"
22import "nx_bigfloat120_div.nx"
23import "nx_bigfloat120_exp.nx"
24import "nx_bigfloat120_ln.nx"
25const K_MAGIC_1100: i64 = 1100
26
27// integer class of y: 0 = non-integer, 1 = odd integer, 2 = even integer
28func pw_int_class(y: i64) -> i64 {
29 let ef: i64 = (y >> 52) & 0x7FF
30 let e: i64 = ef - 1023
31 if e < 0 { return 0 } // |y| < 1: integer only if 0 (caller)
32 if e >= 53 { return 2 } // 2^53+: integer, units bit clear
33 let mant: i64 = (y & 0x000FFFFFFFFFFFFF) | 0x0010000000000000
34 let low: i64 = 52 - e // units bit position
35 if low > 0 {
36 if (mant & ((1 << low) - 1)) != 0 { return 0 }
37 }
38 if ((mant >> low) & 1) == 1 { return 1 }
39 return 2
40}
41
42// |ln(ax)| as bigfloat (out) for finite positive nonzero raw ax; returns the
43// sign of ln (0 = positive, 1 = negative). Structure mirrors bf_ln_f64.
44func bfp_ln_big(out: *i64, ax: i64) -> i64 {
45 var sig: i64 = ax & 0x000FFFFFFFFFFFFF
46 var ef: i64 = (ax >> 52) & 0x7FF
47 if ef == 0 {
48 ef = 1
49 while sig < 0x0010000000000000 { sig = sig << 1; ef = ef - 1 }
50 } else {
51 sig = sig | 0x0010000000000000
52 }
53 var k: i64 = ef - 1023
54 var me: i64 = 0
55 if sig >= bf_sqrt2_cut() { k = k + 1; me = 0 - 1 }
56 let m: *i64 = bf_new()
57 m[0] = me
58 m[1] = sig << 7
59 m[2] = 0
60 let one: *i64 = bf_new()
61 bf_set_int(one, 1)
62 let num: *i64 = bf_new()
63 var sneg: i64 = 0
64 let c: i64 = bf_cmp(m, one)
65 var lnm_zero: i64 = 0
66 if c == 0 {
67 lnm_zero = 1
68 } else {
69 if c > 0 { bf_sub(num, m, one) } else { bf_sub(num, one, m); sneg = 1 }
70 }
71 let lnm: *i64 = bf_new()
72 if lnm_zero == 0 {
73 let den: *i64 = bf_new()
74 bf_add(den, m, one)
75 let s: *i64 = bf_new()
76 bf_div(s, num, den)
77 let z: *i64 = bf_new()
78 bf_mul(z, s, s)
79 let term: *i64 = bf_new()
80 bf_copy(term, s)
81 let sum: *i64 = bf_new()
82 bf_copy(sum, s)
83 let tmul: *i64 = bf_new()
84 let tdiv: *i64 = bf_new()
85 let snew: *i64 = bf_new()
86 var j: i64 = 1
87 while j <= 24 {
88 bf_mul(tmul, term, z)
89 bf_copy(term, tmul)
90 bf_div_small(tdiv, term, 2 * j + 1)
91 bf_add(snew, sum, tdiv)
92 bf_copy(sum, snew)
93 j = j + 1
94 }
95 bf_copy(lnm, sum)
96 lnm[0] = lnm[0] + 1 // * 2: ln(m) = 2*atanh(s)
97 }
98 var ksign: i64 = 0
99 var kabs: i64 = k
100 if k < 0 { ksign = 1; kabs = 0 - k }
101 let kl: *i64 = bf_new()
102 if kabs > 0 {
103 let ln2: *i64 = bf_new()
104 bf_ln2(ln2)
105 bf_mul_small(kl, ln2, kabs)
106 }
107 if lnm_zero == 1 {
108 bf_copy(out, kl) // zero when kabs == 0 (x == 1)
109 return ksign
110 }
111 if kabs == 0 {
112 bf_copy(out, lnm)
113 return sneg
114 }
115 if ksign == sneg {
116 bf_add(out, kl, lnm)
117 return ksign
118 }
119 if bf_cmp(kl, lnm) >= 0 {
120 bf_sub(out, kl, lnm)
121 return ksign
122 }
123 bf_sub(out, lnm, kl)
124 return sneg
125}
126
127// magnitude of exp(+-v) for positive bigfloat v (caller ensures v <= ~1100):
128// k = nearest(v/ln2), r = |v - k*ln2|, 30-term Taylor, 2^k scale; vsign = 1
129// computes the reciprocal. Structure mirrors bf_exp_f64's core.
130func bfp_exp_pm(out: *i64, v: *i64, vsign: i64) -> i64 {
131 let one: *i64 = bf_new()
132 bf_set_int(one, 1)
133 if bf_is_zero(v) == 1 { return bf_copy(out, one) }
134 let ln2: *i64 = bf_new()
135 bf_ln2(ln2)
136 let q: *i64 = bf_new()
137 bf_div(q, v, ln2)
138 let k: i64 = bf_to_int_nearest(q)
139 let kl: *i64 = bf_new()
140 bf_mul_small(kl, ln2, k)
141 let r: *i64 = bf_new()
142 var rneg: i64 = 0
143 if bf_cmp(v, kl) >= 0 {
144 bf_sub(r, v, kl)
145 } else {
146 bf_sub(r, kl, v)
147 rneg = 1
148 }
149 let term: *i64 = bf_new()
150 bf_copy(term, one)
151 let sum: *i64 = bf_new()
152 bf_copy(sum, one)
153 let tmul: *i64 = bf_new()
154 let tdiv: *i64 = bf_new()
155 let snew: *i64 = bf_new()
156 var i: i64 = 1
157 while i <= 30 {
158 bf_mul(tmul, term, r)
159 bf_div_small(tdiv, tmul, i)
160 bf_copy(term, tdiv)
161 bf_add(snew, sum, term)
162 bf_copy(sum, snew)
163 i = i + 1
164 }
165 let er: *i64 = bf_new()
166 if rneg == 0 {
167 bf_copy(er, sum)
168 } else {
169 bf_div(er, one, sum)
170 }
171 er[0] = er[0] + k // * 2^k
172 if vsign == 0 { return bf_copy(out, er) }
173 return bf_div(out, one, er) // exp(-v) = 1/exp(v)
174}
175
176// pow(raw x, raw y) -> raw f64, IEEE-754/C99 semantics.
177func bf_pow_f64(x: i64, y: i64) -> i64 {
178 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
179 let ay: i64 = y & 0x7FFFFFFFFFFFFFFF
180 let xs: i64 = (x >> 63) & 1
181 let ys: i64 = (y >> 63) & 1
182 let one_raw: i64 = 0x3FF0000000000000
183 // pow(x, +-0) = 1 and pow(+1, y) = 1, even for NaN partners
184 if ay == 0 { return one_raw }
185 if x == one_raw { return one_raw }
186 var isnan: i64 = 0
187 if ax > 0x7FF0000000000000 { isnan = 1 }
188 if ay > 0x7FF0000000000000 { isnan = 1 }
189 if isnan == 1 { return 0x7FF8000000000000 }
190 // y = +-inf
191 if ay == 0x7FF0000000000000 {
192 if ax == one_raw { return one_raw } // pow(-1, +-inf) = 1
193 var big: i64 = 0 // result is +inf?
194 if ax > one_raw { big = 1 - ys } else { big = ys }
195 if big == 1 { return 0x7FF0000000000000 }
196 return 0
197 }
198 let ycls: i64 = pw_int_class(ay)
199 // x = +-inf
200 if ax == 0x7FF0000000000000 {
201 var rs: i64 = 0
202 if xs == 1 { if ycls == 1 { rs = 1 } } // -inf with odd-int y keeps sign
203 if ys == 0 { return (rs << 63) | 0x7FF0000000000000 }
204 return rs << 63
205 }
206 // x = +-0
207 if ax == 0 {
208 var rs0: i64 = 0
209 if xs == 1 { if ycls == 1 { rs0 = 1 } }
210 if ys == 0 { return rs0 << 63 } // y > 0: +-0
211 return (rs0 << 63) | 0x7FF0000000000000 // y < 0: +-inf
212 }
213 // negative finite x: integer y only
214 var rsign: i64 = 0
215 if xs == 1 {
216 if ycls == 0 { return 0x7FF8000000000000 }
217 if ycls == 1 { rsign = 1 }
218 }
219 // v = |y * ln|x||, vsign = sign(y*ln|x|)
220 let lt: *i64 = bf_new()
221 let lsign: i64 = bfp_ln_big(lt, ax)
222 if bf_is_zero(lt) == 1 { return (rsign << 63) | one_raw } // |x| = 1
223 let yb: *i64 = bf_new()
224 bf_set_f64(yb, ay)
225 let v: *i64 = bf_new()
226 bf_mul(v, yb, lt)
227 var vsign: i64 = lsign
228 if ys == 1 { vsign = 1 - vsign }
229 // hard band: any |v| > 1100 is decided (finite results end at ~745.2)
230 let band: *i64 = bf_new()
231 bf_set_int(band, K_MAGIC_1100)
232 if bf_cmp(v, band) > 0 {
233 if vsign == 0 { return (rsign << 63) | 0x7FF0000000000000 }
234 return rsign << 63
235 }
236 let mag: *i64 = bf_new()
237 bfp_exp_pm(mag, v, vsign)
238 return bf_to_f64(mag, rsign, 0)
239}