nx_bigfloat120_ln.nx source
↩ module page · 141 lines · 4814 B
1// nx_bigfloat120_ln.nx -- SOVEREIGN ln oracle: binary64-in, binary64-out, 120-bit
2// internals, constants self-computed:
3// sqrt2 boundary = floor(sqrt(2^105)) + 1 via integer digit-pair sqrt (53 pairs,
4// first pair = binary '10', rest zero) -- no imported numbers anywhere.
5// ln(m) = 2*atanh(s), s = (m-1)/(m+1) signed at the boundary; ln(x) = k*ln2 + ln(m).
6// license_tier: ORIGINAL
7
8import "nx_syscalls.nx"
9import "nx_tier.nx"
10import "nx_bigfloat120.nx"
11import "nx_bigfloat120_div.nx"
12import "nx_bigfloat120_exp.nx"
13const K_MAGIC_2047: i64 = 2047
14
15// floor(sqrt(2) * 2^52) + 1: the m-in-[sqrt2/2, sqrt2) split threshold on the
16// 53-bit significand, derived by digit-pair sqrt of the exact radicand 2^105.
17func bf_sqrt2_cut() -> i64 {
18 var rem: i64 = 0
19 var root: i64 = 0
20 var i: i64 = 0
21 while i < 53 {
22 var pair: i64 = 0
23 if i == 0 { pair = 2 } // 2^105 = '10' then 104 zeros
24 rem = (rem << 2) | pair
25 let trial: i64 = (root << 2) | 1
26 if rem >= trial { rem = rem - trial; root = (root << 1) | 1 } else { root = root << 1 }
27 i = i + 1
28 }
29 return root + 1
30}
31
32// ln(raw f64) -> f64 bit pattern.
33func bf_ln_f64(x: i64) -> i64 {
34 let ef0: i64 = (x >> 52) & 0x7FF
35 let sgn: i64 = (x >> 63) & 1
36 if ef0 == K_MAGIC_2047 {
37 if (x & 0x000FFFFFFFFFFFFF) != 0 { return 0x7FF8000000000000 } // NaN
38 if sgn == 1 { return 0x7FF8000000000000 } // ln(-inf)
39 return 0x7FF0000000000000 // ln(+inf)
40 }
41 if (x & 0x7FFFFFFFFFFFFFFF) == 0 { return 0xFFF0000000000000 } // ln(+-0)
42 if sgn == 1 { return 0x7FF8000000000000 } // ln(<0)
43
44 // normalized significand + effective exponent (subnormals pre-normalize)
45 var sig: i64 = x & 0x000FFFFFFFFFFFFF
46 var ef: i64 = ef0
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 // m = sig * 2^(me-52); me in {0,-1}
55 if sig >= bf_sqrt2_cut() { k = k + 1; me = 0 - 1 }
56
57 // m as bigfloat: value = sig * 2^(me-52); with bf e = me, hi = sig<<7
58 let m: *i64 = bf_new()
59 m[0] = me
60 m[1] = sig << 7
61 m[2] = 0
62
63 let one: *i64 = bf_new()
64 bf_set_int(one, 1)
65 // s = (m-1)/(m+1), sign tracked (m < 1 when me = -1 ... or sig = 2^52 exactly)
66 let num: *i64 = bf_new()
67 var sneg: i64 = 0
68 let c: i64 = bf_cmp(m, one)
69 var lnm_zero: i64 = 0
70 if c == 0 {
71 lnm_zero = 1 // ln(1) contribution = 0
72 } else {
73 if c > 0 { bf_sub(num, m, one) } else { bf_sub(num, one, m); sneg = 1 }
74 }
75 let lnm: *i64 = bf_new() // |ln(m)|
76 if lnm_zero == 0 {
77 let den: *i64 = bf_new()
78 bf_add(den, m, one)
79 let s: *i64 = bf_new()
80 bf_div(s, num, den)
81 let z: *i64 = bf_new()
82 bf_mul(z, s, s)
83 let term: *i64 = bf_new()
84 bf_copy(term, s)
85 let sum: *i64 = bf_new()
86 bf_copy(sum, s)
87 let tmul: *i64 = bf_new()
88 let tdiv: *i64 = bf_new()
89 let snew: *i64 = bf_new()
90 var j: i64 = 1
91 while j <= 24 {
92 bf_mul(tmul, term, z)
93 bf_copy(term, tmul)
94 bf_div_small(tdiv, term, 2 * j + 1)
95 bf_add(snew, sum, tdiv)
96 bf_copy(sum, snew)
97 j = j + 1
98 }
99 bf_copy(lnm, sum)
100 lnm[0] = lnm[0] + 1 // * 2: ln(m) = 2*atanh(s)
101 }
102
103 // total = k*ln2 (sign of k) + lnm (sign sneg); combine signed at the boundary
104 var ksign: i64 = 0
105 var kabs: i64 = k
106 if k < 0 { ksign = 1; kabs = 0 - k }
107 let kl: *i64 = bf_new()
108 if kabs > 0 {
109 let ln2: *i64 = bf_new()
110 bf_ln2(ln2)
111 bf_mul_small(kl, ln2, kabs)
112 }
113 // four sign cases collapse: same sign -> add; else subtract smaller from larger
114 var rsign: i64 = 0
115 let tot: *i64 = bf_new()
116 if lnm_zero == 1 {
117 bf_copy(tot, kl)
118 rsign = ksign
119 if kabs == 0 { return 0 } // ln(1) = +0 exactly
120 } else {
121 if kabs == 0 {
122 bf_copy(tot, lnm)
123 rsign = sneg
124 } else {
125 if ksign == sneg {
126 bf_add(tot, kl, lnm)
127 rsign = ksign
128 } else {
129 if bf_cmp(kl, lnm) >= 0 {
130 bf_sub(tot, kl, lnm)
131 rsign = ksign
132 } else {
133 bf_sub(tot, lnm, kl)
134 rsign = sneg
135 }
136 }
137 }
138 }
139 if bf_is_zero(tot) == 1 { return 0 }
140 return bf_to_f64(tot, rsign, 0)
141}