nx_bigfloat120_erfc.nx source
↩ module page · 148 lines · 5161 B
1// nx_bigfloat120_erfc.nx -- ME2-ERF-002 FOUNDATION rung: the SOVEREIGN 120-bit erfc
2// oracle for the 1.75..6 range the erf Taylor series cannot reach (the LOUD-NAN gap the
3// erf rung-1 contract named). No python, no minimax, no external oracle -- and it is
4// PROVABLE by algebra alone (erf(x)+erfc(x)=1), so no trusted reference is imported.
5//
6// Method = the Laplace continued fraction (the classic erfc CF, converges for x>0, fast
7// for x>=1.5):
8// erfc(x) = (exp(-x^2)/sqrt(pi)) * 1/(x + (1/2)/(x + (2/2)/(x + (3/2)/(x + ...))))
9// partial numerators a_k = k/2, partial denominators all = x. Evaluated by backward
10// recurrence from depth N (t=0; for k=N..1: t = a_k/(x+t); CF = 1/(x+t)).
11// exp(-x^2) = 1/exp(x^2); exp(x^2) by DIRECT Taylor (all-positive, x^2<=36 -> the
12// N=160 tail is < 2^-150 relative, no range-reduction needed); sqrt(pi) by Newton.
13// bf_erf here is the SAME 44-term Taylor the trusted erf gate uses, but returns a
14// bigfloat (not f64-rounded) so the gate can check erf+erfc==1 at full precision.
15//
16// All values POSITIVE-magnitude (the bf core convention); erfc and erf are positive on
17// the tested domain so no sign bookkeeping is needed. Composes only proven core ops.
18//
19// module: nishi-core.math.bigfloat120_erfc
20// depends: nishi-core.math.bigfloat120, nishi-core.math.bigfloat120_div, nishi-core.math.bigfloat120_trig
21// capability: SOVEREIGN_ERFC_ORACLE_120BIT
22// license_tier: ORIGINAL
23import "nx_syscalls.nx"
24import "nx_bigfloat120.nx"
25import "nx_bigfloat120_div.nx"
26import "nx_bigfloat120_trig.nx"
27
28// exp(v) for a POSITIVE bigfloat v in [0,36] -- direct Taylor, 160 terms (all positive).
29func bfc_exp(out: *i64, v: *i64) -> i64 {
30 let one: *i64 = bf_new()
31 bf_set_int(one, 1)
32 let term: *i64 = bf_new()
33 bf_copy(term, one)
34 let sum: *i64 = bf_new()
35 bf_copy(sum, one)
36 let tmul: *i64 = bf_new()
37 let tdiv: *i64 = bf_new()
38 let snew: *i64 = bf_new()
39 var i: i64 = 1
40 while i <= 160 {
41 bf_mul(tmul, term, v)
42 bf_div_small(tdiv, tmul, i)
43 bf_copy(term, tdiv)
44 bf_add(snew, sum, term)
45 bf_copy(sum, snew)
46 i = i + 1
47 }
48 bf_copy(out, sum)
49 return 0
50}
51
52// Newton sqrt (8 iterations from an f64 seed; same algorithm the erf/gamma rails use).
53func bfc_sqrt(out: *i64, a: *i64, seed_raw: i64) -> i64 {
54 let s: *i64 = bf_new()
55 bf_set_f64(s, seed_raw)
56 let q: *i64 = bf_new()
57 let sum: *i64 = bf_new()
58 var i: i64 = 0
59 while i < 8 {
60 bf_div(q, a, s)
61 bf_add(sum, s, q)
62 bf_div_small(s, sum, 2)
63 i = i + 1
64 }
65 bf_copy(out, s)
66 return 0
67}
68
69// erfc(x) for a POSITIVE bigfloat xb (x in ~[0.4, 6.5]) via the Laplace CF, depth 240.
70func bf_erfc(out: *i64, xb: *i64) -> i64 {
71 let one: *i64 = bf_new()
72 bf_set_int(one, 1)
73 // x^2 -> exp(x^2) -> exp(-x^2)
74 let x2: *i64 = bf_new()
75 bf_mul(x2, xb, xb)
76 let ex2: *i64 = bf_new()
77 bfc_exp(ex2, x2)
78 let emx2: *i64 = bf_new()
79 bf_div(emx2, one, ex2)
80 // sqrt(pi); prefactor = exp(-x^2)/sqrt(pi)
81 let pi: *i64 = bf_new()
82 bf_pi(pi)
83 let spi: *i64 = bf_new()
84 bfc_sqrt(spi, pi, 0x3FFC5BF891B4EF6A)
85 let pref: *i64 = bf_new()
86 bf_div(pref, emx2, spi)
87 // backward CF: t=0; for k=240..1: t = (k/2)/(x+t)
88 let t: *i64 = bf_new()
89 bf_set_int(t, 0)
90 let ak: *i64 = bf_new()
91 let nk: *i64 = bf_new()
92 let den: *i64 = bf_new()
93 var k: i64 = 240
94 while k >= 1 {
95 bf_set_int(nk, k)
96 bf_div_small(ak, nk, 2) // a_k = k/2
97 bf_add(den, xb, t) // x + t (t = bf zero on the first pass -> x)
98 bf_div(t, ak, den) // t = a_k/(x+t)
99 k = k - 1
100 }
101 bf_add(den, xb, t) // x + t_1
102 let cf: *i64 = bf_new()
103 bf_div(cf, one, den) // 1/(x+t_1)
104 bf_mul(out, pref, cf)
105 return 0
106}
107
108// erf(x) for a POSITIVE bigfloat xb (|x|<=~1.8) via 48-term Taylor, bigfloat output.
109// (the trusted erf-gate oracle; alternating series summed pos/neg then subtracted.)
110func bf_erf(out: *i64, xb: *i64) -> i64 {
111 let pi: *i64 = bf_new()
112 bf_pi(pi)
113 let spi: *i64 = bf_new()
114 bfc_sqrt(spi, pi, 0x3FFC5BF891B4EF6A)
115 let two: *i64 = bf_new()
116 bf_set_int(two, 2)
117 let tosp: *i64 = bf_new()
118 bf_div(tosp, two, spi) // 2/sqrt(pi)
119 let x2: *i64 = bf_new()
120 bf_mul(x2, xb, xb)
121 let u: *i64 = bf_new()
122 bf_copy(u, xb)
123 let pos: *i64 = bf_new()
124 bf_copy(pos, xb) // n=0 term = x (even -> POS)
125 let neg: *i64 = bf_new()
126 bf_set_int(neg, 0)
127 let t: *i64 = bf_new()
128 let term: *i64 = bf_new()
129 let acc: *i64 = bf_new()
130 var n: i64 = 1
131 while n <= 48 {
132 bf_mul(t, u, x2)
133 bf_div_small(u, t, n)
134 bf_div_small(term, u, 2 * n + 1)
135 if (n & 1) == 1 {
136 bf_add(acc, neg, term)
137 bf_copy(neg, acc)
138 } else {
139 bf_add(acc, pos, term)
140 bf_copy(pos, acc)
141 }
142 n = n + 1
143 }
144 let dif: *i64 = bf_new()
145 if bf_cmp(pos, neg) >= 0 { bf_sub(dif, pos, neg) } else { bf_sub(dif, neg, pos) }
146 bf_mul(out, dif, tosp)
147 return 0
148}