nx_bigfloat120_atan.nx source
↩ module page · 141 lines · 5083 B
1// nx_bigfloat120_atan.nx -- SOVEREIGN atan oracle on the 120-bit bigfloat.
2// General atan from NO new primitives: the reduction subtracts table angles
3// atan(1/n) (bf_atan_inv -- the Machin organ already blessed inside bf_pi)
4// via the addition formula
5// atan(x) - atan(1/n) = atan((x - 1/n) / (1 + x/n))
6// choosing n = nearest(1/x) bumped so 1/n <= x: the residual is nonnegative
7// and falls below 1/8 within ~3 rounds from anywhere in [0, 1]; a 30-term
8// Taylor tail finishes (ratio < 2^-6 per term -> truncation < 2^-180).
9// Positive-only core (bigfloat law); signs handled at the f64 boundary.
10// |x| > 1 folds through atan(x) = pi/2 - atan(1/x); atan is total (no refusal).
11//
12// Self-anchor (see _bf_atan_gate_authored): Machin recombination through the
13// GENERAL path reproduces pi, and tan(atan(x)) = x against the blessed
14// bf_sin_r/bf_cos_r organ to < 2^-100 -- no external oracle anywhere.
15// license_tier: ORIGINAL
16
17import "nx_syscalls.nx"
18import "nx_tier.nx"
19import "nx_bigfloat120.nx"
20import "nx_bigfloat120_div.nx"
21import "nx_bigfloat120_trig.nx"
22const K_MAGIC_2047: i64 = 2047
23
24// Taylor tail: out = atan(x) for bigfloat 0 <= x < 1/8 (caller-guaranteed).
25// x - x^3/3 + x^5/5 - ...; pos/neg accumulators (positive-only core).
26func bf_atan_taylor(out: *i64, x: *i64) -> i64 {
27 if bf_is_zero(x) == 1 { return bf_copy(out, x) }
28 let z: *i64 = bf_new()
29 bf_mul(z, x, x)
30 let pw: *i64 = bf_new()
31 bf_copy(pw, x)
32 let pos: *i64 = bf_new()
33 bf_copy(pos, x)
34 let neg: *i64 = bf_new()
35 let tmul: *i64 = bf_new()
36 let tdiv: *i64 = bf_new()
37 let snew: *i64 = bf_new()
38 var k: i64 = 1
39 while k <= 30 {
40 bf_mul(tmul, pw, z)
41 bf_copy(pw, tmul)
42 bf_div_small(tdiv, pw, 2 * k + 1)
43 if (k & 1) == 1 {
44 bf_add(snew, neg, tdiv)
45 bf_copy(neg, snew)
46 } else {
47 bf_add(snew, pos, tdiv)
48 bf_copy(pos, snew)
49 }
50 k = k + 1
51 }
52 return bf_sub(out, pos, neg) // pos > neg (leading term)
53}
54
55// atan on [0, 1]: bigfloat in, bigfloat out. x == 1 -> pi/4 exactly from the
56// Machin pi organ; else subtract table angles until the residual is < 1/8.
57func bf_atan01(out: *i64, xin: *i64) -> i64 {
58 if bf_is_zero(xin) == 1 { return bf_copy(out, xin) }
59 let one: *i64 = bf_new()
60 bf_set_int(one, 1)
61 if bf_cmp(xin, one) == 0 {
62 bf_pi(out)
63 out[0] = out[0] - 2 // pi/4
64 return 0
65 }
66 let eighth: *i64 = bf_new()
67 bf_set_int(eighth, 1)
68 eighth[0] = eighth[0] - 3 // 1/8
69 let x: *i64 = bf_new()
70 bf_copy(x, xin)
71 let acc: *i64 = bf_new()
72 let r: *i64 = bf_new()
73 let invn: *i64 = bf_new()
74 let ang: *i64 = bf_new()
75 let snew: *i64 = bf_new()
76 let num: *i64 = bf_new()
77 let xd: *i64 = bf_new()
78 let den: *i64 = bf_new()
79 let xn: *i64 = bf_new()
80 var round: i64 = 0
81 var go: i64 = 1
82 while go == 1 {
83 if round >= 8 { go = 0 } // mathematical bound is ~3
84 if bf_cmp(x, eighth) <= 0 { go = 0 }
85 if go == 1 {
86 bf_div(r, one, x) // 1/x in (1, 8)
87 var n: i64 = bf_to_int_nearest(r)
88 if n < 2 { n = 2 } // x < 1 here, so n = 1 never valid
89 bf_div_small(invn, one, n)
90 if bf_cmp(invn, x) > 0 { // nearest rounded down: bump once
91 n = n + 1
92 bf_div_small(invn, one, n)
93 }
94 bf_atan_inv(ang, n, 70) // n >= 2: 70 terms > 140 bits
95 bf_add(snew, acc, ang)
96 bf_copy(acc, snew)
97 bf_sub(num, x, invn) // x >= 1/n by construction
98 bf_div_small(xd, x, n)
99 bf_add(den, one, xd)
100 bf_div(xn, num, den)
101 bf_copy(x, xn)
102 round = round + 1
103 }
104 }
105 let tl: *i64 = bf_new()
106 bf_atan_taylor(tl, x)
107 return bf_add(out, acc, tl)
108}
109
110// atan(raw f64) -> f64 bit pattern, round-to-nearest-even. Total domain.
111func bf_atan_f64(x: i64) -> i64 {
112 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
113 let ef: i64 = (x >> 52) & 0x7FF
114 let sgn: i64 = (x >> 63) & 1
115 if ef == K_MAGIC_2047 {
116 if (ax & 0x000FFFFFFFFFFFFF) != 0 { return 0x7FF8000000000000 }
117 let p2: *i64 = bf_new()
118 bf_pi(p2)
119 p2[0] = p2[0] - 1 // atan(+-inf) = +-pi/2
120 return bf_to_f64(p2, sgn, 0)
121 }
122 if ax == 0 { return x } // atan(+-0) = +-0
123 let xb: *i64 = bf_new()
124 bf_set_f64(xb, ax)
125 let one: *i64 = bf_new()
126 bf_set_int(one, 1)
127 let a: *i64 = bf_new()
128 if bf_cmp(xb, one) > 0 {
129 let inv: *i64 = bf_new()
130 bf_div(inv, one, xb)
131 let small: *i64 = bf_new()
132 bf_atan01(small, inv)
133 let p2b: *i64 = bf_new()
134 bf_pi(p2b)
135 p2b[0] = p2b[0] - 1
136 bf_sub(a, p2b, small) // pi/2 > atan01(<1) always
137 return bf_to_f64(a, sgn, 0)
138 }
139 bf_atan01(a, xb)
140 return bf_to_f64(a, sgn, 0)
141}