code wiki / _hdl_build / _bf_gamma_gate_authored.nx
_bf_gamma_gate_authored.nx source
↩ module page · 194 lines · 8338 B
1// _bf_gamma_gate_authored.nx -- self-anchor gate for the SOVEREIGN bigfloat
2// gamma oracle (nx_bigfloat120_gamma). No external oracle:
3// 1. EXACT FACTORIALS: Gamma(n) = (n-1)! for n = 2..19 ((n-1)! <= 18! fits
4// bf_set_int's 2^53 bound), rel < 2^-100 -- Stirling + shifts must land
5// on exact integers across 18 different shift depths.
6// 2. RECURRENCE: Gamma(x+1) = x Gamma(x) at 100 banded non-integer points
7// (x in ~(0.5, 56)): adjacent shift depths must agree, rel < 2^-100.
8// 3. GAMMA(1/2)^2 = pi (squares the result; no sqrt needed), rel < 2^-100.
9// 4. REFLECTION PRODUCT: G(x) G(1-x) sin(pi x) = pi on (0,1) at 40 points
10// (proves the negative-axis path's organs), rel < 2^-95.
11// 5. C99 specials: +-0 -> +-inf, neg integers/-inf -> NaN, +inf -> +inf,
12// 172 -> inf, 171.5 finite, Gamma(1)=Gamma(2)=1 exactly, Gamma(-0.5) vs
13// -2*Gamma(0.5) within 2 ulp, deep-negative signed zeros by parity.
14// Markers: BFGAM-FACT / BFGAM-RECUR / BFGAM-HALF / BFGAM-REFL / BFGAM-SPEC
15
16import "nx_syscalls.nx"
17import "nx_f64.nx"
18import "nx_f64_div.nx"
19import "nx_bigfloat120.nx"
20import "nx_bigfloat120_div.nx"
21import "nx_bigfloat120_trig.nx"
22import "nx_bigfloat120_exp.nx"
23import "nx_bigfloat120_ln.nx"
24import "nx_bigfloat120_pow.nx"
25import "nx_bigfloat120_gamma.nx"
26
27func bgg_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 }
28func bgg_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 }
29func bgg_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 }
30
31func bgg_rng(state: *i64) -> i64 {
32 var s: i64 = state[0]
33 s = s ^ ((s >> 12) & 0x000FFFFFFFFFFFFF)
34 s = s ^ (s << 25)
35 s = s ^ ((s >> 27) & 0x0000001FFFFFFFFF)
36 state[0] = s
37 return s * 2685821657736338717
38}
39
40func bgg_relok(a: *i64, b: *i64, bits: i64) -> i64 {
41 let d: *i64 = bf_new()
42 if bf_cmp(a, b) >= 0 { bf_sub(d, a, b) } else { bf_sub(d, b, a) }
43 if bf_is_zero(d) == 1 { return 1 }
44 if d[0] <= b[0] - bits { return 1 }
45 return 0
46}
47
48func bgg_ord(v: i64) -> i64 {
49 if v >= 0 { return v }
50 return (1 << 63) - v
51}
52
53func bgg_ulp(a: i64, b: i64) -> i64 {
54 if a == b { return 0 }
55 var u: i64 = bgg_ord(a) - bgg_ord(b)
56 if u < 0 { u = 0 - u }
57 return u
58}
59
60func main() -> i64 {
61 var bad: i64 = 0
62 let one: *i64 = bf_new()
63 bf_set_int(one, 1)
64 let st: *i64 = sys_mmap(16) as *i64
65 st[0] = 17724538509055160273
66
67 // ---- 1. exact factorials across shift depths ----
68 var fbad: i64 = 0
69 var fact: i64 = 1
70 var n: i64 = 2
71 while n <= 19 {
72 fact = fact * (n - 1)
73 let xb: *i64 = bf_new()
74 bf_set_int(xb, n)
75 let g: *i64 = bf_new()
76 bfg_gamma_pos(g, xb)
77 let want: *i64 = bf_new()
78 bf_set_int(want, fact)
79 if bgg_relok(g, want, 100) == 0 {
80 fbad = fbad + 1
81 bgg_puts("BFGAM-FACTBAD n=" as *u8); bgg_putn(n); bgg_puts("\n" as *u8)
82 }
83 n = n + 1
84 }
85 bgg_puts("BFGAM-FACT total=18 bad=" as *u8); bgg_putn(fbad); bgg_puts(" bar=rel<2^-100\n" as *u8)
86 bad = bad + fbad
87
88 // ---- 2. recurrence at banded non-integer points ----
89 var rbad: i64 = 0
90 var j: i64 = 0
91 while j < 100 {
92 var raw: i64 = bgg_rng(st)
93 let efr: i64 = 1021 + (raw & 7) // 0.25 .. 64
94 raw = (efr << 52) | (bgg_rng(st) & 0x000FFFFFFFFFFFFF)
95 let xb: *i64 = bf_new()
96 bf_set_f64(xb, raw)
97 let ga: *i64 = bf_new()
98 bfg_gamma_pos(ga, xb)
99 let xb1: *i64 = bf_new()
100 bf_add(xb1, xb, one)
101 let gb: *i64 = bf_new()
102 bfg_gamma_pos(gb, xb1)
103 let xa: *i64 = bf_new()
104 bf_mul(xa, xb, ga)
105 if bgg_relok(xa, gb, 100) == 0 {
106 rbad = rbad + 1
107 if rbad <= 10 { bgg_puts("BFGAM-RECURBAD x=" as *u8); bgg_puthex(raw); bgg_puts("\n" as *u8) }
108 }
109 j = j + 1
110 }
111 bgg_puts("BFGAM-RECUR total=100 bad=" as *u8); bgg_putn(rbad); bgg_puts(" bar=rel<2^-100\n" as *u8)
112 bad = bad + rbad
113
114 // ---- 3. Gamma(1/2)^2 = pi ----
115 let halfb: *i64 = bf_new()
116 bf_copy(halfb, one)
117 halfb[0] = halfb[0] - 1
118 let gh: *i64 = bf_new()
119 bfg_gamma_pos(gh, halfb)
120 let gh2: *i64 = bf_new()
121 bf_mul(gh2, gh, gh)
122 let pi: *i64 = bf_new()
123 bf_pi(pi)
124 if bgg_relok(gh2, pi, 100) == 1 {
125 bgg_puts("BFGAM-HALF ok Gamma(1/2)^2=pi rel<2^-100\n" as *u8)
126 } else {
127 bad = bad + 1
128 bgg_puts("BFGAM-HALF BAD\n" as *u8)
129 }
130
131 // ---- 4. reflection product on (0,1) ----
132 var xbad: i64 = 0
133 j = 0
134 while j < 40 {
135 var raw2: i64 = bgg_rng(st)
136 let efr2: i64 = 1015 + (raw2 & 7) // 2^-8 .. 1
137 raw2 = (efr2 << 52) | (bgg_rng(st) & 0x000FFFFFFFFFFFFF)
138 let xb2: *i64 = bf_new()
139 bf_set_f64(xb2, raw2)
140 let g1: *i64 = bf_new()
141 bfg_gamma_pos(g1, xb2)
142 let omx: *i64 = bf_new()
143 bf_sub(omx, one, xb2) // x < 1 here
144 let g2: *i64 = bf_new()
145 bfg_gamma_pos(g2, omx)
146 let px: *i64 = bf_new()
147 bf_mul(px, pi, xb2)
148 let sv: *i64 = bf_new()
149 bfg_sin_big(sv, px) // sin(pi x) > 0 on (0,1)
150 let p1: *i64 = bf_new()
151 bf_mul(p1, g1, g2)
152 let p2: *i64 = bf_new()
153 bf_mul(p2, p1, sv)
154 if bgg_relok(p2, pi, 95) == 0 {
155 xbad = xbad + 1
156 if xbad <= 5 { bgg_puts("BFGAM-REFLBAD x=" as *u8); bgg_puthex(raw2); bgg_puts("\n" as *u8) }
157 }
158 j = j + 1
159 }
160 bgg_puts("BFGAM-REFL total=40 bad=" as *u8); bgg_putn(xbad); bgg_puts(" bar=rel<2^-95\n" as *u8)
161 bad = bad + xbad
162
163 // ---- 5. C99 specials ----
164 var sbad: i64 = 0
165 if bf_gamma_f64(0) != 0x7FF0000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD a\n" as *u8) }
166 if bf_gamma_f64(1 << 63) != 0xFFF0000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD b\n" as *u8) }
167 if bf_gamma_f64(0x7FF8000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD c\n" as *u8) }
168 if bf_gamma_f64(0x7FF0000000000000) != 0x7FF0000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD d\n" as *u8) }
169 if bf_gamma_f64(0xFFF0000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD e\n" as *u8) }
170 if bf_gamma_f64(0xC008000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD f\n" as *u8) } // -3
171 if bf_gamma_f64(0x4065800000000000) != 0x7FF0000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD g\n" as *u8) } // 172
172 let fin: i64 = bf_gamma_f64(0x40656FFFFFFFFFFF) // ~171.49: finite
173 if ((fin >> 52) & 0x7FF) == 2047 { sbad = sbad + 1; bgg_puts("SPECBAD h\n" as *u8) }
174 if bf_gamma_f64(0x3FF0000000000000) != 0x3FF0000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD i\n" as *u8) } // G(1)=1
175 if bf_gamma_f64(0x4000000000000000) != 0x3FF0000000000000 { sbad = sbad + 1; bgg_puts("SPECBAD j\n" as *u8) } // G(2)=1
176 // Gamma(-1/2) = -2 Gamma(1/2): exact x2, compare within 2 ulp
177 let gph: i64 = bf_gamma_f64(0x3FE0000000000000)
178 let gmh: i64 = bf_gamma_f64(0xBFE0000000000000)
179 let want2: i64 = (gph + (1 << 52)) | (1 << 63) // -2*Gamma(1/2)
180 if bgg_ulp(gmh, want2) > 2 { sbad = sbad + 1; bgg_puts("SPECBAD k\n" as *u8) }
181 // deep-negative signed zeros by interval parity
182 if bf_gamma_f64(0xC069080000000000) != (1 << 63) { sbad = sbad + 1; bgg_puts("SPECBAD l\n" as *u8) } // -200.25: n even -> -0
183 if bf_gamma_f64(0xC069280000000000) != 0 { sbad = sbad + 1; bgg_puts("SPECBAD m\n" as *u8) } // -201.25: n odd -> +0
184 bgg_puts("BFGAM-SPEC bad=" as *u8); bgg_putn(sbad); bgg_puts("\n" as *u8)
185 bad = bad + sbad
186
187 if bad == 0 {
188 bgg_puts("BFGAM-GATE verdict=GREEN\n" as *u8)
189 return 0
190 }
191 bgg_puts("BFGAM-GATE verdict=RED\n" as *u8)
192 if bad > 100 { bad = 100 }
193 return bad
194}