code wiki / _hdl_build / _bf_trig_gate_authored.nx
_bf_trig_gate_authored.nx source
↩ module page · 134 lines · 5425 B
1// _bf_trig_gate_authored.nx -- gate for the SOVEREIGN trig oracle. All anchors
2// are internal (no python, no external constants):
3// P0 pi by TWO independent formulas agrees to ~2^-117:
4// Machin: pi = 16*atan(1/5) - 4*atan(1/239)
5// Gauss: pi = 48*atan(1/18) + 32*atan(1/57) - 20*atan(1/239)
6// P1 Pythagorean identity in bigfloat at 200 xorshift points (|x| < 2^20):
7// |sin^2(r)+cos^2(r) - 1| < 2^-110 after full reduction
8// P2 f64-level IEEE invariants: sin odd / cos even bit-exact, sin(+-0)=+-0,
9// cos(+-0)=1, NaN/inf -> NaN, |x| >= 2^20 -> NaN (honest v1 refusal)
10// Markers: TRIGP0 ... / TRIGP1 ... / TRIGP2 ... / BF-TRIG-GATE verdict=
11
12import "nx_syscalls.nx"
13import "nx_bigfloat120.nx"
14import "nx_bigfloat120_div.nx"
15import "nx_bigfloat120_trig.nx"
16
17func bft_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
18func bft_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 }
19
20func bft_rng(state: *i64) -> i64 {
21 var s: i64 = state[0]
22 s = s ^ ((s >> 12) & 0x000FFFFFFFFFFFFF)
23 s = s ^ (s << 25)
24 s = s ^ ((s >> 27) & 0x0000001FFFFFFFFF)
25 state[0] = s
26 return s * 2685821657736338717
27}
28
29func main() -> i64 {
30 var bad: i64 = 0
31
32 // ---- P0: cross-formula pi ----
33 let pi1: *i64 = bf_new()
34 bf_pi(pi1)
35 let a18: *i64 = bf_new()
36 bf_atan_inv(a18, 18, 16)
37 let a57: *i64 = bf_new()
38 bf_atan_inv(a57, 57, 12)
39 let a239: *i64 = bf_new()
40 bf_atan_inv(a239, 239, 9)
41 let t1: *i64 = bf_new()
42 bf_mul_small(t1, a18, 48)
43 let t2: *i64 = bf_new()
44 bf_mul_small(t2, a57, 32)
45 let t3: *i64 = bf_new()
46 bf_mul_small(t3, a239, 20)
47 let pi2: *i64 = bf_new()
48 bf_add(pi2, t1, t2)
49 let pi2b: *i64 = bf_new()
50 bf_sub(pi2b, pi2, t3)
51 let dif: *i64 = bf_new()
52 if bf_cmp(pi1, pi2b) >= 0 { bf_sub(dif, pi1, pi2b) } else { bf_sub(dif, pi2b, pi1) }
53 var p0ok: i64 = 0
54 if bf_is_zero(dif) == 1 { p0ok = 1 }
55 if p0ok == 0 {
56 if dif[0] <= pi1[0] - 115 { p0ok = 1 }
57 }
58 if p0ok == 0 { bad = bad + 1 }
59 bft_puts("TRIGP0 machin-vs-gauss agree=" as *u8); bft_putn(p0ok)
60 if p0ok == 0 { bft_puts(" diff-e=" as *u8); bft_putn(dif[0]) }
61 bft_puts("\n" as *u8)
62
63 // ---- P1: pythagorean identity in-bigfloat ----
64 let one: *i64 = bf_new()
65 bf_set_int(one, 1)
66 let st: *i64 = sys_mmap(16) as *i64
67 st[0] = 424242424242421
68 var p1bad: i64 = 0
69 var i: i64 = 0
70 while i < 200 {
71 // random f64 in |x| < 2^20: exponent field 1..1042
72 var raw: i64 = bft_rng(st)
73 var efr: i64 = 800 + (raw & 255) // 2^-223 .. 2^-... wide band
74 if (raw & 1024) != 0 { efr = 1023 + (raw & 15) } // and a near-1..2^16 band
75 raw = (efr << 52) | (bft_rng(st) & 0x000FFFFFFFFFFFFF)
76 let xb: *i64 = bf_new()
77 bf_set_f64(xb, raw)
78 let r: *i64 = bf_new()
79 let kq: *i64 = sys_mmap(16) as *i64
80 _bf_trig_reduce(r, kq, xb)
81 let sv: *i64 = bf_new()
82 let cv: *i64 = bf_new()
83 bf_sin_r(sv, r)
84 bf_cos_r(cv, r)
85 let s2: *i64 = bf_new()
86 bf_mul(s2, sv, sv)
87 let c2: *i64 = bf_new()
88 bf_mul(c2, cv, cv)
89 let sum: *i64 = bf_new()
90 bf_add(sum, s2, c2)
91 let d2: *i64 = bf_new()
92 if bf_cmp(sum, one) >= 0 { bf_sub(d2, sum, one) } else { bf_sub(d2, one, sum) }
93 var ok: i64 = 0
94 if bf_is_zero(d2) == 1 { ok = 1 }
95 if ok == 0 { if d2[0] <= 0 - 110 { ok = 1 } }
96 if ok == 0 { p1bad = p1bad + 1 }
97 i = i + 1
98 }
99 if p1bad > 0 { bad = bad + 1 }
100 bft_puts("TRIGP1 identity points=200 bad=" as *u8); bft_putn(p1bad); bft_puts("\n" as *u8)
101
102 // ---- P2: f64 IEEE invariants ----
103 var p2bad: i64 = 0
104 st[0] = 777000111222333
105 var j: i64 = 0
106 while j < 150 {
107 var raw2: i64 = bft_rng(st)
108 var ef2: i64 = 950 + (raw2 & 63)
109 raw2 = (ef2 << 52) | (bft_rng(st) & 0x000FFFFFFFFFFFFF)
110 let sp: i64 = bf_sin_f64(raw2)
111 let sn: i64 = bf_sin_f64(raw2 | (1 << 63))
112 if sn != (sp ^ (1 << 63)) { p2bad = p2bad + 1 } // sin odd
113 let cp: i64 = bf_cos_f64(raw2)
114 let cn: i64 = bf_cos_f64(raw2 | (1 << 63))
115 if cn != cp { p2bad = p2bad + 1 } // cos even
116 j = j + 1
117 }
118 if bf_sin_f64(0) != 0 { p2bad = p2bad + 1 }
119 if bf_sin_f64(1 << 63) != (1 << 63) { p2bad = p2bad + 1 }
120 if bf_cos_f64(0) != 0x3FF0000000000000 { p2bad = p2bad + 1 }
121 if bf_cos_f64(1 << 63) != 0x3FF0000000000000 { p2bad = p2bad + 1 }
122 if bf_sin_f64(0x7FF0000000000000) != 0x7FF8000000000000 { p2bad = p2bad + 1 }
123 if bf_cos_f64(0x7FF8000000000000) != 0x7FF8000000000000 { p2bad = p2bad + 1 }
124 if bf_sin_f64(0x4140000000000000) != 0x7FF8000000000000 { p2bad = p2bad + 1 } // 2^21 refused
125 if p2bad > 0 { bad = bad + 1 }
126 bft_puts("TRIGP2 invariants bad=" as *u8); bft_putn(p2bad); bft_puts("\n" as *u8)
127
128 if bad == 0 {
129 bft_puts("BF-TRIG-GATE verdict=GREEN\n" as *u8)
130 return 0
131 }
132 bft_puts("BF-TRIG-GATE verdict=RED\n" as *u8)
133 return bad
134}