code wiki / _hdl_build / _bf_sinhcosh_gate_authored.nx
_bf_sinhcosh_gate_authored.nx source
↩ module page · 191 lines · 8469 B
1// _bf_sinhcosh_gate_authored.nx -- self-anchor gate for the SOVEREIGN bigfloat
2// sinh/cosh oracle (nx_bigfloat120_sinh). No external oracle:
3// 1. DOUBLE-ANGLE (cancellation-free, full domain): sinh(2u) = 2 sinh u cosh u
4// and cosh(2u) = 1 + 2 sinh^2 u at 150 banded points (u < 350), rel < 2^-100.
5// 2. PYTHAGOREAN strip: cosh^2 - sinh^2 = 1 for |x| <= 5 (the subtraction
6// cancels 2*log2(cosh) <= 13 bits there), 60 points, rel < 2^-90.
7// 3. f64 tie to the BLESSED exp organ: sinh(x)+cosh(x) within 2 ulp of
8// bf_exp_f64(x) at 100 points (different exp code path).
9// 4. Specials + overflow boundary + Taylor/exp seam at x = 1.
10// Markers: BFSH-DOUBLE / BFSH-PYTH / BFSH-EXPTIE / BFSH-SPEC / verdict=
11
12import "nx_syscalls.nx"
13import "nx_f64.nx"
14import "nx_f64_div.nx"
15import "nx_bigfloat120.nx"
16import "nx_bigfloat120_div.nx"
17import "nx_bigfloat120_exp.nx"
18import "nx_bigfloat120_ln.nx"
19import "nx_bigfloat120_pow.nx"
20import "nx_bigfloat120_sinh.nx"
21
22func bsg_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
23func bsg_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 }
24func bsg_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 }
25
26func bsg_rng(state: *i64) -> i64 {
27 var s: i64 = state[0]
28 s = s ^ ((s >> 12) & 0x000FFFFFFFFFFFFF)
29 s = s ^ (s << 25)
30 s = s ^ ((s >> 27) & 0x0000001FFFFFFFFF)
31 state[0] = s
32 return s * 2685821657736338717
33}
34
35func bsg_relok(a: *i64, b: *i64, bits: i64) -> i64 {
36 let d: *i64 = bf_new()
37 if bf_cmp(a, b) >= 0 { bf_sub(d, a, b) } else { bf_sub(d, b, a) }
38 if bf_is_zero(d) == 1 { return 1 }
39 if d[0] <= b[0] - bits { return 1 }
40 return 0
41}
42
43func bsg_ord(v: i64) -> i64 {
44 if v >= 0 { return v }
45 return (1 << 63) - v
46}
47
48func bsg_ulp(a: i64, b: i64) -> i64 {
49 if a == b { return 0 }
50 var u: i64 = bsg_ord(a) - bsg_ord(b)
51 if u < 0 { u = 0 - u }
52 return u
53}
54
55func main() -> i64 {
56 var bad: i64 = 0
57 let one: *i64 = bf_new()
58 bf_set_int(one, 1)
59 let st: *i64 = sys_mmap(16) as *i64
60 st[0] = 31622776601683793320
61
62 // ---- 1. double-angle, full domain ----
63 var dbad: i64 = 0
64 var j: i64 = 0
65 while j < 150 {
66 var raw: i64 = bsg_rng(st)
67 var efr: i64 = 0
68 let band: i64 = raw & 3
69 if band == 0 { efr = 1015 + (raw & 15) } // 2^-8 .. 2^8
70 if band == 1 { efr = 1023 + (raw & 7) } // 1 .. 256
71 if band == 2 { efr = 990 + (raw & 31) } // small
72 if band == 3 { efr = 1029 + (raw & 3) } // 64 .. 512 (2u < 710)
73 raw = (efr << 52) | (bsg_rng(st) & 0x000FFFFFFFFFFFFF)
74 if efr == 1031 { raw = (1031 << 52) | (bsg_rng(st) & 0x0007FFFFFFFFFFFF) } // cap < 350
75 if efr == 1032 { raw = (1029 << 52) | (bsg_rng(st) & 0x000FFFFFFFFFFFFF) }
76 let u: *i64 = bf_new()
77 bf_set_f64(u, raw)
78 let sh: *i64 = bf_new()
79 bf_sinh_xb(sh, u)
80 let ch: *i64 = bf_new()
81 bf_cosh_xb(ch, u)
82 let u2: *i64 = bf_new()
83 bf_copy(u2, u)
84 u2[0] = u2[0] + 1 // exact 2u
85 let sh2: *i64 = bf_new()
86 bf_sinh_xb(sh2, u2)
87 let ch2: *i64 = bf_new()
88 bf_cosh_xb(ch2, u2)
89 // sinh(2u) = 2 sinh u cosh u
90 let pr: *i64 = bf_new()
91 bf_mul(pr, sh, ch)
92 pr[0] = pr[0] + 1
93 var ok: i64 = bsg_relok(pr, sh2, 100)
94 // cosh(2u) = 1 + 2 sinh^2 u
95 let s2: *i64 = bf_new()
96 bf_mul(s2, sh, sh)
97 s2[0] = s2[0] + 1
98 let rhs: *i64 = bf_new()
99 bf_add(rhs, one, s2)
100 if bsg_relok(rhs, ch2, 100) == 0 { ok = 0 }
101 if ok == 0 {
102 dbad = dbad + 1
103 if dbad <= 10 { bsg_puts("BFSH-DBLBAD u=" as *u8); bsg_puthex(raw); bsg_puts("\n" as *u8) }
104 }
105 j = j + 1
106 }
107 bsg_puts("BFSH-DOUBLE total=150 bad=" as *u8); bsg_putn(dbad); bsg_puts(" bar=rel<2^-100\n" as *u8)
108 bad = bad + dbad
109
110 // ---- 2. pythagorean strip |x| <= 5 ----
111 var ybad: i64 = 0
112 j = 0
113 while j < 60 {
114 var raw2: i64 = bsg_rng(st)
115 let efr2: i64 = 1018 + (raw2 & 7) // 2^-5 .. 8 capped
116 raw2 = (efr2 << 52) | (bsg_rng(st) & 0x000FFFFFFFFFFFFF)
117 if efr2 == 1025 { raw2 = (1024 << 52) | (bsg_rng(st) & 0x000FFFFFFFFFFFFF) } // < 4
118 let xv: *i64 = bf_new()
119 bf_set_f64(xv, raw2)
120 let sh3: *i64 = bf_new()
121 bf_sinh_xb(sh3, xv)
122 let ch3: *i64 = bf_new()
123 bf_cosh_xb(ch3, xv)
124 let cq: *i64 = bf_new()
125 bf_mul(cq, ch3, ch3)
126 let sq: *i64 = bf_new()
127 bf_mul(sq, sh3, sh3)
128 let df: *i64 = bf_new()
129 bf_sub(df, cq, sq) // cosh^2 > sinh^2
130 if bsg_relok(df, one, 90) == 0 {
131 ybad = ybad + 1
132 if ybad <= 5 { bsg_puts("BFSH-PYTHBAD x=" as *u8); bsg_puthex(raw2); bsg_puts("\n" as *u8) }
133 }
134 j = j + 1
135 }
136 bsg_puts("BFSH-PYTH total=60 bad=" as *u8); bsg_putn(ybad); bsg_puts(" bar=rel<2^-90\n" as *u8)
137 bad = bad + ybad
138
139 // ---- 3. f64 tie: sinh + cosh = e^x vs the blessed exp organ ----
140 var ebad: i64 = 0
141 j = 0
142 while j < 100 {
143 var raw3: i64 = bsg_rng(st)
144 let efr3: i64 = 1015 + (raw3 & 15)
145 raw3 = (efr3 << 52) | (bsg_rng(st) & 0x000FFFFFFFFFFFFF)
146 let sf: i64 = bf_sinh_f64(raw3)
147 let cf: i64 = bf_cosh_f64(raw3)
148 let ex: i64 = bf_exp_f64(raw3)
149 if bsg_ulp(nx_f64_add(sf, cf), ex) > 2 {
150 ebad = ebad + 1
151 if ebad <= 5 { bsg_puts("BFSH-EXPBAD x=" as *u8); bsg_puthex(raw3); bsg_puts("\n" as *u8) }
152 }
153 j = j + 1
154 }
155 bsg_puts("BFSH-EXPTIE total=100 bad=" as *u8); bsg_putn(ebad); bsg_puts(" bar=2ulp\n" as *u8)
156 bad = bad + ebad
157
158 // ---- 4. specials + overflow boundary + seam ----
159 var sbad: i64 = 0
160 if bf_sinh_f64(0) != 0 { sbad = sbad + 1; bsg_puts("SPECBAD a\n" as *u8) }
161 if bf_sinh_f64(1 << 63) != (1 << 63) { sbad = sbad + 1; bsg_puts("SPECBAD b\n" as *u8) }
162 if bf_cosh_f64(0) != 0x3FF0000000000000 { sbad = sbad + 1; bsg_puts("SPECBAD c\n" as *u8) }
163 if bf_cosh_f64(1 << 63) != 0x3FF0000000000000 { sbad = sbad + 1; bsg_puts("SPECBAD d\n" as *u8) }
164 if bf_sinh_f64(0x7FF8000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; bsg_puts("SPECBAD e\n" as *u8) }
165 if bf_sinh_f64(0x7FF0000000000000) != 0x7FF0000000000000 { sbad = sbad + 1; bsg_puts("SPECBAD f\n" as *u8) }
166 if bf_sinh_f64(0xFFF0000000000000) != 0xFFF0000000000000 { sbad = sbad + 1; bsg_puts("SPECBAD g\n" as *u8) }
167 if bf_cosh_f64(0xFFF0000000000000) != 0x7FF0000000000000 { sbad = sbad + 1; bsg_puts("SPECBAD h\n" as *u8) }
168 // overflow boundary: true cut ~710.476: 710.0 finite, 710.5 -> inf
169 if bf_sinh_f64(0x4086340000000000) != 0x7FF0000000000000 { sbad = sbad + 1; bsg_puts("SPECBAD i\n" as *u8) } // 710.5 -> inf
170 let fb: i64 = bf_sinh_f64(0x4086300000000000) // 710.0: finite
171 if ((fb >> 52) & 0x7FF) == 2047 { sbad = sbad + 1; bsg_puts("SPECBAD i2\n" as *u8) }
172 let fin: i64 = bf_cosh_f64(0x4086200000000000) // 708.25: finite
173 if ((fin >> 52) & 0x7FF) == 2047 { sbad = sbad + 1; bsg_puts("SPECBAD j\n" as *u8) }
174 if bf_cosh_f64(0x40C0000000000000) != 0x7FF0000000000000 { sbad = sbad + 1; bsg_puts("SPECBAD k\n" as *u8) } // 8192 -> inf
175 // sinh odd at the Taylor/exp seam
176 let pv: i64 = bf_sinh_f64(0x3FF0000000000000)
177 if bf_sinh_f64(0xBFF0000000000000) != (pv | (1 << 63)) { sbad = sbad + 1; bsg_puts("SPECBAD l\n" as *u8) }
178 // monotone across the seam x = 1
179 let s0: i64 = bf_sinh_f64(0x3FEFFFFFFFFFFFFF)
180 if s0 > pv { sbad = sbad + 1; bsg_puts("SPECBAD m\n" as *u8) }
181 bsg_puts("BFSH-SPEC bad=" as *u8); bsg_putn(sbad); bsg_puts("\n" as *u8)
182 bad = bad + sbad
183
184 if bad == 0 {
185 bsg_puts("BFSH-GATE verdict=GREEN\n" as *u8)
186 return 0
187 }
188 bsg_puts("BFSH-GATE verdict=RED\n" as *u8)
189 if bad > 100 { bad = 100 }
190 return bad
191}