nx_bigfloat120_exp.nx source
↩ module page · 125 lines · 4103 B
1// nx_bigfloat120_exp.nx -- SOVEREIGN exp oracle: binary64-in, binary64-out, all
2// internals 120-bit bigfloat, EVERY constant self-computed (no mpmath, no python):
3// ln2 = 2 * atanh(1/3) = 2 * (1/3 + (1/3)^3/3 + (1/3)^5/5 + ...) [39 odd terms]
4// Range reduction k = nearest(x / ln2) is taken in bigfloat; r = |x - k*ln2| with
5// the sign handled by the reciprocal trick (exp(-r) = 1/exp(r)) so the Taylor
6// series sums positive terms only. 30 terms: r <= ln2/2 gives a 2^-122 tail.
7// license_tier: ORIGINAL
8
9import "nx_syscalls.nx"
10import "nx_tier.nx"
11import "nx_bigfloat120.nx"
12import "nx_bigfloat120_div.nx"
13const K_MAGIC_2047: i64 = 2047
14const K_MAGIC_1033: i64 = 1033
15
16// ln2 to 120 bits, computed on demand (callers cache the block).
17func bf_ln2(out: *i64) -> i64 {
18 let third: *i64 = bf_new()
19 let one: *i64 = bf_new()
20 bf_set_int(one, 1)
21 bf_div_small(third, one, 3) // 1/3
22 let z: *i64 = bf_new()
23 bf_mul(z, third, third) // (1/3)^2 = 1/9
24 let term: *i64 = bf_new()
25 bf_copy(term, third)
26 let sum: *i64 = bf_new()
27 bf_copy(sum, third)
28 let tnext: *i64 = bf_new()
29 let tdiv: *i64 = bf_new()
30 let snew: *i64 = bf_new()
31 var j: i64 = 1
32 while j <= 39 {
33 bf_mul(tnext, term, z)
34 bf_copy(term, tnext)
35 bf_div_small(tdiv, term, 2 * j + 1)
36 bf_add(snew, sum, tdiv)
37 bf_copy(sum, snew)
38 j = j + 1
39 }
40 bf_copy(out, sum)
41 out[0] = out[0] + 1 // * 2 (atanh -> ln2)
42 return 0
43}
44
45// exp(raw f64) -> f64 bit pattern, sovereign 120-bit pipeline.
46func bf_exp_f64(x: i64) -> i64 {
47 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
48 let ef: i64 = (x >> 52) & 0x7FF
49 let sgn: i64 = (x >> 63) & 1
50 // specials
51 if ef == K_MAGIC_2047 {
52 if (x & 0x000FFFFFFFFFFFFF) != 0 { return 0x7FF8000000000000 } // NaN
53 if sgn == 1 { return 0 } // exp(-inf)
54 return 0x7FF0000000000000 // exp(+inf)
55 }
56 if ax == 0 { return 0x3FF0000000000000 } // exp(0)=1
57 // range cuts (|x| > 746 covers both overflow and deep underflow)
58 if ef >= K_MAGIC_1033 {
59 // |x| >= 1024: way past both cuts
60 if sgn == 1 { return 0 }
61 return 0x7FF0000000000000
62 }
63 let ln2: *i64 = bf_new()
64 bf_ln2(ln2)
65 let xb: *i64 = bf_new()
66 bf_set_f64(xb, ax)
67 // overflow/underflow precise cuts via bigfloat compare: 710 and 746
68 let cut: *i64 = bf_new()
69 if sgn == 0 {
70 bf_set_int(cut, 710)
71 if bf_cmp(xb, cut) > 0 { return 0x7FF0000000000000 }
72 } else {
73 bf_set_int(cut, 746)
74 if bf_cmp(xb, cut) > 0 { return 0 }
75 }
76 // k = nearest(|x| / ln2), then r = | |x| - k*ln2 |
77 let q: *i64 = bf_new()
78 bf_div(q, xb, ln2)
79 let k: i64 = bf_to_int_nearest(q)
80 let kl: *i64 = bf_new()
81 bf_mul_small(kl, ln2, k)
82 let r: *i64 = bf_new()
83 var rneg: i64 = 0 // r sign within |x| frame
84 if bf_cmp(xb, kl) >= 0 {
85 bf_sub(r, xb, kl)
86 } else {
87 bf_sub(r, kl, xb)
88 rneg = 1
89 }
90 // Taylor exp(r), positive r, 30 terms
91 let one: *i64 = bf_new()
92 bf_set_int(one, 1)
93 let term: *i64 = bf_new()
94 bf_copy(term, one)
95 let sum: *i64 = bf_new()
96 bf_copy(sum, one)
97 let tmul: *i64 = bf_new()
98 let tdiv: *i64 = bf_new()
99 let snew: *i64 = bf_new()
100 var i: i64 = 1
101 while i <= 30 {
102 bf_mul(tmul, term, r)
103 bf_div_small(tdiv, tmul, i)
104 bf_copy(term, tdiv)
105 bf_add(snew, sum, term)
106 bf_copy(sum, snew)
107 i = i + 1
108 }
109 // exp(|x|) = 2^k * sum (rneg=0) or 2^k / sum (rneg=1)
110 let er: *i64 = bf_new()
111 if rneg == 0 {
112 bf_copy(er, sum)
113 } else {
114 bf_div(er, one, sum)
115 }
116 var kk: i64 = k
117 // exp(x) for negative x = 1 / exp(|x|) = 2^-k / er'
118 if sgn == 1 {
119 let inv: *i64 = bf_new()
120 bf_div(inv, one, er)
121 bf_copy(er, inv)
122 kk = 0 - k
123 }
124 return bf_to_f64(er, 0, kk)
125}