nx_bigfloat120_trig.nx source
↩ module page · 218 lines · 7556 B
1// nx_bigfloat120_trig.nx -- SOVEREIGN trig oracle on the 120-bit bigfloat: pi by
2// Machin's formula, sin/cos by pi/2 reduction + Taylor. No imported constants:
3// pi = 16*atan(1/5) - 4*atan(1/239) (Machin 1706; series self-computed)
4// Alternating series run on TWO positive accumulators (core is positive-only),
5// subtracted once at the end -- no catastrophic cancellation (neg/pos ratio is
6// bounded well below 1 for 1/n arguments and r <= pi/4).
7//
8// Self-anchor: the trig gate's phase 0 checks sin^2 + cos^2 = 1 to < 2^-110 at
9// random points -- an internal-consistency proof that needs no external oracle.
10//
11// v1 domain: |x| < 2^20 (the f64 kernel's 3-part Cody-Waite validity range);
12// huge-argument reduction (Payne-Hanek class) is a named follow-on rung.
13// license_tier: ORIGINAL
14
15import "nx_syscalls.nx"
16import "nx_tier.nx"
17import "nx_bigfloat120.nx"
18import "nx_bigfloat120_div.nx"
19const K_MAGIC_57121: i64 = 57121
20const K_MAGIC_2047: i64 = 2047
21const K_MAGIC_1043: i64 = 1043
22
23// atan(1/n) for small integer n: sum_k (-1)^k * (1/n)^(2k+1) / (2k+1).
24// Positive and negative terms accumulate separately.
25func bf_atan_inv(out: *i64, n: i64, terms: i64) -> i64 {
26 let one: *i64 = bf_new()
27 bf_set_int(one, 1)
28 let x: *i64 = bf_new()
29 bf_div_small(x, one, n) // 1/n
30 let z: *i64 = bf_new()
31 bf_mul(z, x, x) // 1/n^2
32 let term: *i64 = bf_new()
33 bf_copy(term, x)
34 let pos: *i64 = bf_new()
35 bf_copy(pos, x) // k = 0 term is positive
36 let neg: *i64 = bf_new()
37 let tmul: *i64 = bf_new()
38 let tdiv: *i64 = bf_new()
39 let snew: *i64 = bf_new()
40 var k: i64 = 1
41 while k <= terms {
42 bf_mul(tmul, term, z)
43 bf_copy(term, tmul)
44 bf_div_small(tdiv, term, 2 * k + 1)
45 if (k & 1) == 1 {
46 bf_add(snew, neg, tdiv)
47 bf_copy(neg, snew)
48 } else {
49 bf_add(snew, pos, tdiv)
50 bf_copy(pos, snew)
51 }
52 k = k + 1
53 }
54 return bf_sub(out, pos, neg) // pos > neg always (leading term)
55}
56
57// pi to 120 bits: 16*atan(1/5) - 4*atan(1/239)
58func bf_pi(out: *i64) -> i64 {
59 let a5: *i64 = bf_new()
60 bf_atan_inv(a5, 5, 28) // ratio 1/25: 28 terms ~ 2^-134
61 let a239: *i64 = bf_new()
62 bf_atan_inv(a239, 239, 9) // ratio 1/K_MAGIC_57121: 9 terms ~ 2^-150
63 a5[0] = a5[0] + 4 // * 16
64 a239[0] = a239[0] + 2 // * 4
65 return bf_sub(out, a5, a239)
66}
67
68// sin(r) on bigfloat, 0 <= r <= pi/4: r - r^3/3! + r^5/5! - ... (17+17 terms)
69func bf_sin_r(out: *i64, r: *i64) -> i64 {
70 if bf_is_zero(r) == 1 { return bf_copy(out, r) }
71 let z: *i64 = bf_new()
72 bf_mul(z, r, r)
73 let term: *i64 = bf_new()
74 bf_copy(term, r)
75 let pos: *i64 = bf_new()
76 bf_copy(pos, r)
77 let neg: *i64 = bf_new()
78 let tmul: *i64 = bf_new()
79 let tdiv: *i64 = bf_new()
80 let snew: *i64 = bf_new()
81 var k: i64 = 1
82 while k <= 16 {
83 bf_mul(tmul, term, z)
84 bf_div_small(tdiv, tmul, (2 * k) * (2 * k + 1))
85 bf_copy(term, tdiv)
86 if (k & 1) == 1 {
87 bf_add(snew, neg, term)
88 bf_copy(neg, snew)
89 } else {
90 bf_add(snew, pos, term)
91 bf_copy(pos, snew)
92 }
93 k = k + 1
94 }
95 return bf_sub(out, pos, neg)
96}
97
98// cos(r) on bigfloat, 0 <= r <= pi/4: 1 - r^2/2! + r^4/4! - ...
99func bf_cos_r(out: *i64, r: *i64) -> i64 {
100 let one: *i64 = bf_new()
101 bf_set_int(one, 1)
102 if bf_is_zero(r) == 1 { return bf_copy(out, one) }
103 let z: *i64 = bf_new()
104 bf_mul(z, r, r)
105 let term: *i64 = bf_new()
106 bf_copy(term, one)
107 let pos: *i64 = bf_new()
108 bf_copy(pos, one)
109 let neg: *i64 = bf_new()
110 let tmul: *i64 = bf_new()
111 let tdiv: *i64 = bf_new()
112 let snew: *i64 = bf_new()
113 var k: i64 = 1
114 while k <= 16 {
115 bf_mul(tmul, term, z)
116 bf_div_small(tdiv, tmul, (2 * k - 1) * (2 * k))
117 bf_copy(term, tdiv)
118 if (k & 1) == 1 {
119 bf_add(snew, neg, term)
120 bf_copy(neg, snew)
121 } else {
122 bf_add(snew, pos, term)
123 bf_copy(pos, snew)
124 }
125 k = k + 1
126 }
127 return bf_sub(out, pos, neg)
128}
129
130// shared reduction: |x| -> quadrant k (mod 4) + r in [0, pi/4-ish] with the
131// half-quadrant fold. Writes k&3 to kq[0], |r| to r, r-sign flag to kq[1]
132// (sin(q*pio2 + s*r) expansion handles the fold via the quadrant tables below).
133// Approach: k = nearest(|x| / (pi/2)); rr = |x| - k*(pi/2) signed.
134func _bf_trig_reduce(r: *i64, kq: *i64, xabs: *i64) -> i64 {
135 let pi: *i64 = bf_new()
136 bf_pi(pi)
137 let pio2: *i64 = bf_new()
138 bf_copy(pio2, pi)
139 pio2[0] = pio2[0] - 1
140 let q: *i64 = bf_new()
141 bf_div(q, xabs, pio2)
142 let k: i64 = bf_to_int_nearest(q)
143 let kl: *i64 = bf_new()
144 bf_mul_small(kl, pio2, k)
145 var rneg: i64 = 0
146 if bf_cmp(xabs, kl) >= 0 {
147 bf_sub(r, xabs, kl)
148 } else {
149 bf_sub(r, kl, xabs)
150 rneg = 1
151 }
152 kq[0] = k & 3
153 kq[1] = rneg
154 return 0
155}
156
157// sin(raw f64) -> f64 bit pattern (|x| < 2^20; NaN past that = honest refusal)
158func bf_sin_f64(x: i64) -> i64 {
159 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
160 let ef: i64 = (x >> 52) & 0x7FF
161 let sgn: i64 = (x >> 63) & 1
162 if ef == K_MAGIC_2047 { return 0x7FF8000000000000 } // inf and NaN -> NaN
163 if ax == 0 { return x } // sin(+-0) = +-0
164 if ef >= K_MAGIC_1043 { return 0x7FF8000000000000 } // |x| >= 2^20: refused v1
165 let xb: *i64 = bf_new()
166 bf_set_f64(xb, ax)
167 let r: *i64 = bf_new()
168 let kq: *i64 = sys_mmap(16) as *i64
169 _bf_trig_reduce(r, kq, xb)
170 let sv: *i64 = bf_new()
171 let cv: *i64 = bf_new()
172 bf_sin_r(sv, r)
173 bf_cos_r(cv, r)
174 // sin(k*pio2 + t) where t = +-r: quadrant table
175 // k=0: +-sin(r) k=1: cos(r) [t-sign irrelevant for magnitude via
176 // k=2: -+sin(r) cos(-r)=cos(r)] k=3: -cos(r)
177 var outsig: i64 = 0
178 var usecos: i64 = 0
179 let kk: i64 = kq[0]
180 let rneg: i64 = kq[1]
181 if kk == 0 { outsig = rneg }
182 if kk == 1 { usecos = 1; outsig = 0 }
183 if kk == 2 { outsig = 1 - rneg }
184 if kk == 3 { usecos = 1; outsig = 1 }
185 if sgn == 1 { outsig = 1 - outsig } // sin odd
186 if usecos == 1 { return bf_to_f64(cv, outsig, 0) }
187 return bf_to_f64(sv, outsig, 0)
188}
189
190// cos(raw f64) -> f64 bit pattern (|x| < 2^20)
191func bf_cos_f64(x: i64) -> i64 {
192 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
193 let ef: i64 = (x >> 52) & 0x7FF
194 if ef == K_MAGIC_2047 { return 0x7FF8000000000000 }
195 if ax == 0 { return 0x3FF0000000000000 } // cos(+-0) = 1
196 if ef >= K_MAGIC_1043 { return 0x7FF8000000000000 }
197 let xb: *i64 = bf_new()
198 bf_set_f64(xb, ax)
199 let r: *i64 = bf_new()
200 let kq: *i64 = sys_mmap(16) as *i64
201 _bf_trig_reduce(r, kq, xb)
202 let sv: *i64 = bf_new()
203 let cv: *i64 = bf_new()
204 bf_sin_r(sv, r)
205 bf_cos_r(cv, r)
206 // cos(k*pio2 + t), t = +-r: k=0: cos k=1: -+sin (-sin(t), t=+-r)
207 // k=2: -cos k=3: +-sin
208 var outsig: i64 = 0
209 var usesin: i64 = 0
210 let kk: i64 = kq[0]
211 let rneg: i64 = kq[1]
212 if kk == 0 { outsig = 0 }
213 if kk == 1 { usesin = 1; outsig = 1 - rneg }
214 if kk == 2 { outsig = 1 }
215 if kk == 3 { usesin = 1; outsig = rneg }
216 if usesin == 1 { return bf_to_f64(sv, outsig, 0) } // cos even: x sign irrelevant
217 return bf_to_f64(cv, outsig, 0)
218}