code wiki / _hdl_build / _bf_tan_gate_authored.nx
_bf_tan_gate_authored.nx source
↩ module page · 195 lines · 8337 B
1// _bf_tan_gate_authored.nx -- self-anchor gate for the SOVEREIGN bigfloat tan
2// oracle (nx_bigfloat120_tan). No external oracle:
3// 1. HALF-ANGLE: tan(2u) = 2t/(1-t^2), t = tan(u), at 150 banded points --
4// u and 2u take DIFFERENT reduction paths (k vs ~2k), so the algebraic
5// identity catches reduction and quadrant-table defects. rel < 2^-100
6// (points with |1-t^2| below 2^-30 are skipped-and-counted: cancellation
7// there honestly exceeds the bar).
8// 2. PI-SHIFT: tan(x + pi) = tan(x) with the EXACT bigfloat pi, 50 points.
9// 3. f64 tie to the BLESSED organs: |tan(x) - fl(sin(x)/cos(x))| <= 2 ulp
10// (quotient of two correctly-rounded values bounds at ~1.5 ulp + 0.5).
11// 4. Specials + pole signs: tan(fl(pi/2)) huge POSITIVE, next ulp NEGATIVE,
12// oddness, +-0, NaN/inf/refusal.
13// Markers: BFTAN-HALF / BFTAN-PISHIFT / BFTAN-QUOT / BFTAN-SPEC / verdict=
14
15import "nx_syscalls.nx"
16import "nx_f64.nx"
17import "nx_f64_div.nx"
18import "nx_bigfloat120.nx"
19import "nx_bigfloat120_div.nx"
20import "nx_bigfloat120_trig.nx"
21import "nx_bigfloat120_tan.nx"
22
23func btg_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
24func btg_putn(v: i64) -> i64 { let bb: *u8 = sys_mmap(28); var m: i64 = v; if m < 0 { m = 0 - m; sys_write(1, "-" as *u8, 1) }; let t: *u8 = sys_mmap(28); var k: i64 = 0; if m == 0 { t[0] = 48; k = 1 }; while m > 0 { t[k] = 48 + (m % 10); m = m / 10; k = k + 1 }; var i: i64 = 0; while i < k { bb[i] = t[k-1-i]; i = i + 1 }; sys_write(1, bb, k); return 0 }
25func btg_puthex(v: i64) -> i64 { let bb: *u8 = sys_mmap(20); var i: i64 = 0; while i < 16 { let nib: i64 = (v >> ((15 - i) * 4)) & 15; if nib < 10 { bb[i] = 48 + nib } else { bb[i] = 55 + nib } i = i + 1 } sys_write(1, bb, 16); return 0 }
26
27func btg_rng(state: *i64) -> i64 {
28 var s: i64 = state[0]
29 s = s ^ ((s >> 12) & 0x000FFFFFFFFFFFFF)
30 s = s ^ (s << 25)
31 s = s ^ ((s >> 27) & 0x0000001FFFFFFFFF)
32 state[0] = s
33 return s * 2685821657736338717
34}
35
36func btg_relok(a: *i64, b: *i64, bits: i64) -> i64 {
37 let d: *i64 = bf_new()
38 if bf_cmp(a, b) >= 0 { bf_sub(d, a, b) } else { bf_sub(d, b, a) }
39 if bf_is_zero(d) == 1 { return 1 }
40 if d[0] <= b[0] - bits { return 1 }
41 return 0
42}
43
44func btg_ord(v: i64) -> i64 {
45 if v >= 0 { return v }
46 return (1 << 63) - v
47}
48
49func btg_ulp(a: i64, b: i64) -> i64 {
50 if a == b { return 0 }
51 var u: i64 = btg_ord(a) - btg_ord(b)
52 if u < 0 { u = 0 - u }
53 return u
54}
55
56func main() -> i64 {
57 var bad: i64 = 0
58 let one: *i64 = bf_new()
59 bf_set_int(one, 1)
60 let st: *i64 = sys_mmap(16) as *i64
61 st[0] = 26457513110645905905
62
63 // ---- 1. half-angle across reduction paths ----
64 var hbad: i64 = 0
65 var skipped: i64 = 0
66 var j: i64 = 0
67 while j < 150 {
68 var raw: i64 = btg_rng(st)
69 var efr: i64 = 0
70 let band: i64 = raw & 3
71 if band == 0 { efr = 1021 + (raw & 3) } // 0.25 .. 4
72 if band == 1 { efr = 1023 + (raw & 15) } // 1 .. 2^16
73 if band == 2 { efr = 1015 + (raw & 7) } // small
74 if band == 3 { efr = 1033 + (raw & 7) } // 2^10 .. 2^18
75 raw = (efr << 52) | (btg_rng(st) & 0x000FFFFFFFFFFFFF)
76 let xb: *i64 = bf_new()
77 bf_set_f64(xb, raw)
78 let m1: *i64 = bf_new()
79 let s1: i64 = bf_tan_xb(m1, xb)
80 let xb2: *i64 = bf_new()
81 bf_copy(xb2, xb)
82 xb2[0] = xb2[0] + 1 // exact 2u
83 let m2: *i64 = bf_new()
84 let s2: i64 = bf_tan_xb(m2, xb2)
85 // rhs = 2*m1 / |1 - m1^2|, sign s1 ^ (m1 > 1)
86 let msq: *i64 = bf_new()
87 bf_mul(msq, m1, m1)
88 let den: *i64 = bf_new()
89 var dsg: i64 = 0
90 if bf_cmp(msq, one) >= 0 { bf_sub(den, msq, one); dsg = 1 } else { bf_sub(den, one, msq) }
91 var skip: i64 = 0
92 if bf_is_zero(den) == 1 { skip = 1 }
93 if skip == 0 { if den[0] < 0 - 30 { skip = 1 } }
94 if skip == 1 {
95 skipped = skipped + 1
96 } else {
97 let num: *i64 = bf_new()
98 bf_copy(num, m1)
99 num[0] = num[0] + 1 // 2*m1
100 let rhs: *i64 = bf_new()
101 bf_div(rhs, num, den)
102 // reduction carries ~2^(e_u - 117) ABSOLUTE error; when r lands
103 // small that dominates -- scale the bar by u's exponent.
104 var bits: i64 = 100
105 if xb2[0] > 0 { bits = 100 - xb2[0] }
106 if bits < 60 { bits = 60 }
107 var ok: i64 = btg_relok(rhs, m2, bits)
108 if (s1 ^ dsg) != s2 { ok = 0 }
109 if ok == 0 {
110 hbad = hbad + 1
111 if hbad <= 10 {
112 btg_puts("BFTAN-HALFBAD u=" as *u8); btg_puthex(raw); btg_puts("\n" as *u8)
113 }
114 }
115 }
116 j = j + 1
117 }
118 btg_puts("BFTAN-HALF total=150 bad=" as *u8); btg_putn(hbad)
119 btg_puts(" skipped=" as *u8); btg_putn(skipped)
120 btg_puts(" bar=rel<2^-100\n" as *u8)
121 bad = bad + hbad
122
123 // ---- 2. pi-shift invariance with the exact bigfloat pi ----
124 let pi: *i64 = bf_new()
125 bf_pi(pi)
126 var pbad: i64 = 0
127 j = 0
128 while j < 50 {
129 var raw2: i64 = btg_rng(st)
130 let efr2: i64 = 1020 + (raw2 & 15)
131 raw2 = (efr2 << 52) | (btg_rng(st) & 0x000FFFFFFFFFFFFF)
132 let xa: *i64 = bf_new()
133 bf_set_f64(xa, raw2)
134 let ma: *i64 = bf_new()
135 let sa: i64 = bf_tan_xb(ma, xa)
136 let xc: *i64 = bf_new()
137 bf_add(xc, xa, pi)
138 let mc: *i64 = bf_new()
139 let sc: i64 = bf_tan_xb(mc, xc)
140 var pok: i64 = btg_relok(mc, ma, 100)
141 if sa != sc { pok = 0 }
142 if pok == 0 {
143 pbad = pbad + 1
144 if pbad <= 5 { btg_puts("BFTAN-PIBAD x=" as *u8); btg_puthex(raw2); btg_puts("\n" as *u8) }
145 }
146 j = j + 1
147 }
148 btg_puts("BFTAN-PISHIFT total=50 bad=" as *u8); btg_putn(pbad); btg_puts("\n" as *u8)
149 bad = bad + pbad
150
151 // ---- 3. f64 tie to blessed sin/cos: quotient within 2 ulp ----
152 var qbad: i64 = 0
153 j = 0
154 while j < 100 {
155 var raw3: i64 = btg_rng(st)
156 let efr3: i64 = 1018 + (raw3 & 15)
157 raw3 = (efr3 << 52) | (btg_rng(st) & 0x000FFFFFFFFFFFFF)
158 if (btg_rng(st) & 1) == 1 { raw3 = raw3 | (1 << 63) }
159 let tg: i64 = bf_tan_f64(raw3)
160 let qq: i64 = nx_f64_div(bf_sin_f64(raw3), bf_cos_f64(raw3))
161 if btg_ulp(tg, qq) > 2 {
162 qbad = qbad + 1
163 if qbad <= 5 { btg_puts("BFTAN-QUOTBAD x=" as *u8); btg_puthex(raw3); btg_puts("\n" as *u8) }
164 }
165 j = j + 1
166 }
167 btg_puts("BFTAN-QUOT total=100 bad=" as *u8); btg_putn(qbad); btg_puts(" bar=2ulp-of-blessed-quotient\n" as *u8)
168 bad = bad + qbad
169
170 // ---- 4. specials + pole signs ----
171 var sbad: i64 = 0
172 if bf_tan_f64(0) != 0 { sbad = sbad + 1; btg_puts("SPECBAD a\n" as *u8) }
173 if bf_tan_f64(1 << 63) != (1 << 63) { sbad = sbad + 1; btg_puts("SPECBAD b\n" as *u8) }
174 if bf_tan_f64(0x7FF8000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; btg_puts("SPECBAD c\n" as *u8) }
175 if bf_tan_f64(0x7FF0000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; btg_puts("SPECBAD d\n" as *u8) }
176 if bf_tan_f64(0x4140000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; btg_puts("SPECBAD e\n" as *u8) } // 2^21 refused
177 let below: i64 = bf_tan_f64(0x3FF921FB54442D18) // fl(pi/2) < pi/2: huge +
178 let above: i64 = bf_tan_f64(0x3FF921FB54442D19) // next ulp: crossed pole, -
179 if ((below >> 63) & 1) != 0 { sbad = sbad + 1; btg_puts("SPECBAD f\n" as *u8) }
180 if ((below >> 52) & 0x7FF) < 1073 { sbad = sbad + 1; btg_puts("SPECBAD g\n" as *u8) } // |tan| > 2^50
181 if ((above >> 63) & 1) != 1 { sbad = sbad + 1; btg_puts("SPECBAD h\n" as *u8) }
182 if ((above >> 52) & 0x7FF) < 1073 { sbad = sbad + 1; btg_puts("SPECBAD i\n" as *u8) }
183 let pv: i64 = bf_tan_f64(0x3FE0000000000000)
184 if bf_tan_f64(0xBFE0000000000000) != (pv | (1 << 63)) { sbad = sbad + 1; btg_puts("SPECBAD j\n" as *u8) }
185 btg_puts("BFTAN-SPEC bad=" as *u8); btg_putn(sbad); btg_puts("\n" as *u8)
186 bad = bad + sbad
187
188 if bad == 0 {
189 btg_puts("BFTAN-GATE verdict=GREEN\n" as *u8)
190 return 0
191 }
192 btg_puts("BFTAN-GATE verdict=RED\n" as *u8)
193 if bad > 100 { bad = 100 }
194 return bad
195}