code wiki / _hdl_build / _f64_sinhcosh_gate_authored.nx
_f64_sinhcosh_gate_authored.nx source
↩ module page · 155 lines · 6252 B
1// _f64_sinhcosh_gate_authored.nx -- ME1 sinh/cosh gate: dd kernels
2// (_pe_f64sinhcosh) vs the LIVE sovereign bigfloat oracle. Pure Nishi.
3// Edges: Taylor/exp seam at 1, I-drop seam at 28, overflow boundary
4// (710.0 / ln(2^1025) / 710.5), tiny/subnormal; bands across (0,1), [1,28),
5// [28,710]. GREEN iff max ULP <= 1 on BOTH functions.
6// Markers: SHBAD/CHBAD x= exp= got= ulp= / SINH-ULP / COSH-ULP / verdict=
7
8import "nx_syscalls.nx"
9import "nx_f64.nx"
10import "nx_f64_div.nx"
11import "nx_f64_cvt.nx"
12import "nx_bigfloat120.nx"
13import "nx_bigfloat120_div.nx"
14import "nx_bigfloat120_exp.nx"
15import "nx_bigfloat120_ln.nx"
16import "nx_bigfloat120_pow.nx"
17import "nx_bigfloat120_sinh.nx"
18import "_pm_sinhcosh_spec.nx"
19import "_pe_f64sinhcosh.nx"
20
21func shg_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
22func shg_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 }
23func shg_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 }
24
25func shg_rng(state: *i64) -> i64 {
26 var s: i64 = state[0]
27 s = s ^ ((s >> 12) & 0x000FFFFFFFFFFFFF)
28 s = s ^ (s << 25)
29 s = s ^ ((s >> 27) & 0x0000001FFFFFFFFF)
30 state[0] = s
31 return s * 2685821657736338717
32}
33
34func shg_ord(v: i64) -> i64 {
35 if v >= 0 { return v }
36 return (1 << 63) - v
37}
38
39func shg_ulp(got: i64, want: i64) -> i64 {
40 if got == want { return 0 }
41 var u: i64 = shg_ord(got) - shg_ord(want)
42 if u < 0 { u = 0 - u }
43 return u
44}
45
46func main() -> i64 {
47 let xs: *i64 = sys_mmap(8 * 512) as *i64
48 var n: i64 = 0
49 xs[n] = 0x7FF8000000000000; n = n + 1
50 xs[n] = 0x7FF0000000000000; n = n + 1
51 xs[n] = 0xFFF0000000000000; n = n + 1
52 xs[n] = 0; n = n + 1
53 xs[n] = 1 << 63; n = n + 1
54 xs[n] = 1; n = n + 1 // min subnormal
55 xs[n] = 0x000FFFFFFFFFFFFF; n = n + 1
56 xs[n] = 0x3E50000000000000; n = n + 1 // 2^-26
57 xs[n] = 0x3FEFFFFFFFFFFFFF; n = n + 1 // 1 - ulp (Taylor side)
58 xs[n] = 0x3FF0000000000000; n = n + 1 // 1.0 (exp side)
59 xs[n] = 0x3FF0000000000001; n = n + 1
60 xs[n] = 0x403BFFFFFFFFFFFF; n = n + 1 // 28 - ulp (I-branch)
61 xs[n] = 0x403C000000000000; n = n + 1 // 28 (E/2 branch)
62 xs[n] = 0x4086300000000000; n = n + 1 // 710.0 finite
63 xs[n] = 0x40862E42FEFA39EF; n = n + 1 // ln(2^1025): the brink
64 xs[n] = 0x4086340000000000; n = n + 1 // 710.5 inf
65 xs[n] = 0x4086380000000000; n = n + 1 // 711 inf
66 xs[n] = 0xC086300000000000; n = n + 1 // -710 (sinh sign)
67 xs[n] = 0xBFE0000000000000; n = n + 1 // -0.5
68 let st: *i64 = sys_mmap(16) as *i64
69 st[0] = 27182818284590452353
70 var j: i64 = 0
71 while j < 320 {
72 var raw: i64 = shg_rng(st)
73 var efr: i64 = 0
74 let band: i64 = raw & 3
75 if band == 0 { efr = 1015 + (raw & 7) } // (2^-8, 2) Taylor zone
76 if band == 1 { efr = 1023 + (raw & 3) } // [1, 16) I-branch
77 if band == 2 { efr = 1027 + (raw & 3) } // [16, 256) seam + E/2
78 if band == 3 { efr = 990 + (raw & 63) } // tiny through mid
79 var rw: i64 = (efr << 52) | (shg_rng(st) & 0x000FFFFFFFFFFFFF)
80 if efr == 1032 { rw = (1032 << 52) | (shg_rng(st) & 0x0006FFFFFFFFFFFF) } // cap < ~700
81 if (shg_rng(st) & 1) == 1 { rw = rw | (1 << 63) }
82 xs[n] = rw
83 n = n + 1
84 j = j + 1
85 }
86
87 var s_exact: i64 = 0
88 var s_u1: i64 = 0
89 var s_worse: i64 = 0
90 var s_max: i64 = 0
91 var c_exact: i64 = 0
92 var c_u1: i64 = 0
93 var c_worse: i64 = 0
94 var c_max: i64 = 0
95 var i: i64 = 0
96 while i < n {
97 let x: i64 = xs[i]
98 let sw: i64 = bf_sinh_f64(x)
99 let sg: i64 = nx_f64_sinh(x)
100 let su: i64 = shg_ulp(sg, sw)
101 if su == 0 { s_exact = s_exact + 1 }
102 if su == 1 { s_u1 = s_u1 + 1 }
103 if su > 1 {
104 s_worse = s_worse + 1
105 if s_worse <= 20 {
106 shg_puts("SHBAD x=" as *u8); shg_puthex(x)
107 shg_puts(" exp=" as *u8); shg_puthex(sw)
108 shg_puts(" got=" as *u8); shg_puthex(sg)
109 shg_puts(" ulp=" as *u8); shg_putn(su)
110 shg_puts("\n" as *u8)
111 }
112 }
113 if su > s_max { s_max = su }
114 let cw: i64 = bf_cosh_f64(x)
115 let cg: i64 = nx_f64_cosh(x)
116 let cu: i64 = shg_ulp(cg, cw)
117 if cu == 0 { c_exact = c_exact + 1 }
118 if cu == 1 { c_u1 = c_u1 + 1 }
119 if cu > 1 {
120 c_worse = c_worse + 1
121 if c_worse <= 20 {
122 shg_puts("CHBAD x=" as *u8); shg_puthex(x)
123 shg_puts(" exp=" as *u8); shg_puthex(cw)
124 shg_puts(" got=" as *u8); shg_puthex(cg)
125 shg_puts(" ulp=" as *u8); shg_putn(cu)
126 shg_puts("\n" as *u8)
127 }
128 }
129 if cu > c_max { c_max = cu }
130 i = i + 1
131 }
132
133 shg_puts("SINH-ULP total=" as *u8); shg_putn(n)
134 shg_puts(" exact=" as *u8); shg_putn(s_exact)
135 shg_puts(" ulp1=" as *u8); shg_putn(s_u1)
136 shg_puts(" worse=" as *u8); shg_putn(s_worse)
137 shg_puts(" max=" as *u8); shg_putn(s_max)
138 shg_puts("\n" as *u8)
139 shg_puts("COSH-ULP total=" as *u8); shg_putn(n)
140 shg_puts(" exact=" as *u8); shg_putn(c_exact)
141 shg_puts(" ulp1=" as *u8); shg_putn(c_u1)
142 shg_puts(" worse=" as *u8); shg_putn(c_worse)
143 shg_puts(" max=" as *u8); shg_putn(c_max)
144 shg_puts("\n" as *u8)
145 if s_worse == 0 {
146 if c_worse == 0 {
147 shg_puts("SINHCOSH-GATE verdict=GREEN\n" as *u8)
148 return 0
149 }
150 }
151 shg_puts("SINHCOSH-GATE verdict=RED\n" as *u8)
152 var rc: i64 = s_worse + c_worse
153 if rc > 100 { rc = 100 }
154 return rc
155}