code wiki / _hdl_build / _bf_tan_gate_authored.nx

_bf_tan_gate_authored.nx source

↩ module page · 198 lines · 8853 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 // LN40 (2026-09-03): this PRNG seed was written as a 20-digit decimal that does not fit 64 bits, so 62 // every build before today WRAPPED it silently. The value below is exactly what it wrapped to, so the 63 // sequence this gate has always exercised is unchanged -- the literal now says what the machine does. 64 st[0] = 8010769036936354289 65 66 // ---- 1. half-angle across reduction paths ---- 67 var hbad: i64 = 0 68 var skipped: i64 = 0 69 var j: i64 = 0 70 while j < 150 { 71 var raw: i64 = btg_rng(st) 72 var efr: i64 = 0 73 let band: i64 = raw & 3 74 if band == 0 { efr = 1021 + (raw & 3) } // 0.25 .. 4 75 if band == 1 { efr = 1023 + (raw & 15) } // 1 .. 2^16 76 if band == 2 { efr = 1015 + (raw & 7) } // small 77 if band == 3 { efr = 1033 + (raw & 7) } // 2^10 .. 2^18 78 raw = (efr << 52) | (btg_rng(st) & 0x000FFFFFFFFFFFFF) 79 let xb: *i64 = bf_new() 80 bf_set_f64(xb, raw) 81 let m1: *i64 = bf_new() 82 let s1: i64 = bf_tan_xb(m1, xb) 83 let xb2: *i64 = bf_new() 84 bf_copy(xb2, xb) 85 xb2[0] = xb2[0] + 1 // exact 2u 86 let m2: *i64 = bf_new() 87 let s2: i64 = bf_tan_xb(m2, xb2) 88 // rhs = 2*m1 / |1 - m1^2|, sign s1 ^ (m1 > 1) 89 let msq: *i64 = bf_new() 90 bf_mul(msq, m1, m1) 91 let den: *i64 = bf_new() 92 var dsg: i64 = 0 93 if bf_cmp(msq, one) >= 0 { bf_sub(den, msq, one); dsg = 1 } else { bf_sub(den, one, msq) } 94 var skip: i64 = 0 95 if bf_is_zero(den) == 1 { skip = 1 } 96 if skip == 0 { if den[0] < 0 - 30 { skip = 1 } } 97 if skip == 1 { 98 skipped = skipped + 1 99 } else { 100 let num: *i64 = bf_new() 101 bf_copy(num, m1) 102 num[0] = num[0] + 1 // 2*m1 103 let rhs: *i64 = bf_new() 104 bf_div(rhs, num, den) 105 // reduction carries ~2^(e_u - 117) ABSOLUTE error; when r lands 106 // small that dominates -- scale the bar by u's exponent. 107 var bits: i64 = 100 108 if xb2[0] > 0 { bits = 100 - xb2[0] } 109 if bits < 60 { bits = 60 } 110 var ok: i64 = btg_relok(rhs, m2, bits) 111 if (s1 ^ dsg) != s2 { ok = 0 } 112 if ok == 0 { 113 hbad = hbad + 1 114 if hbad <= 10 { 115 btg_puts("BFTAN-HALFBAD u=" as *u8); btg_puthex(raw); btg_puts("\n" as *u8) 116 } 117 } 118 } 119 j = j + 1 120 } 121 btg_puts("BFTAN-HALF total=150 bad=" as *u8); btg_putn(hbad) 122 btg_puts(" skipped=" as *u8); btg_putn(skipped) 123 btg_puts(" bar=rel<2^-100\n" as *u8) 124 bad = bad + hbad 125 126 // ---- 2. pi-shift invariance with the exact bigfloat pi ---- 127 let pi: *i64 = bf_new() 128 bf_pi(pi) 129 var pbad: i64 = 0 130 j = 0 131 while j < 50 { 132 var raw2: i64 = btg_rng(st) 133 let efr2: i64 = 1020 + (raw2 & 15) 134 raw2 = (efr2 << 52) | (btg_rng(st) & 0x000FFFFFFFFFFFFF) 135 let xa: *i64 = bf_new() 136 bf_set_f64(xa, raw2) 137 let ma: *i64 = bf_new() 138 let sa: i64 = bf_tan_xb(ma, xa) 139 let xc: *i64 = bf_new() 140 bf_add(xc, xa, pi) 141 let mc: *i64 = bf_new() 142 let sc: i64 = bf_tan_xb(mc, xc) 143 var pok: i64 = btg_relok(mc, ma, 100) 144 if sa != sc { pok = 0 } 145 if pok == 0 { 146 pbad = pbad + 1 147 if pbad <= 5 { btg_puts("BFTAN-PIBAD x=" as *u8); btg_puthex(raw2); btg_puts("\n" as *u8) } 148 } 149 j = j + 1 150 } 151 btg_puts("BFTAN-PISHIFT total=50 bad=" as *u8); btg_putn(pbad); btg_puts("\n" as *u8) 152 bad = bad + pbad 153 154 // ---- 3. f64 tie to blessed sin/cos: quotient within 2 ulp ---- 155 var qbad: i64 = 0 156 j = 0 157 while j < 100 { 158 var raw3: i64 = btg_rng(st) 159 let efr3: i64 = 1018 + (raw3 & 15) 160 raw3 = (efr3 << 52) | (btg_rng(st) & 0x000FFFFFFFFFFFFF) 161 if (btg_rng(st) & 1) == 1 { raw3 = raw3 | (1 << 63) } 162 let tg: i64 = bf_tan_f64(raw3) 163 let qq: i64 = nx_f64_div(bf_sin_f64(raw3), bf_cos_f64(raw3)) 164 if btg_ulp(tg, qq) > 2 { 165 qbad = qbad + 1 166 if qbad <= 5 { btg_puts("BFTAN-QUOTBAD x=" as *u8); btg_puthex(raw3); btg_puts("\n" as *u8) } 167 } 168 j = j + 1 169 } 170 btg_puts("BFTAN-QUOT total=100 bad=" as *u8); btg_putn(qbad); btg_puts(" bar=2ulp-of-blessed-quotient\n" as *u8) 171 bad = bad + qbad 172 173 // ---- 4. specials + pole signs ---- 174 var sbad: i64 = 0 175 if bf_tan_f64(0) != 0 { sbad = sbad + 1; btg_puts("SPECBAD a\n" as *u8) } 176 if bf_tan_f64(1 << 63) != (1 << 63) { sbad = sbad + 1; btg_puts("SPECBAD b\n" as *u8) } 177 if bf_tan_f64(0x7FF8000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; btg_puts("SPECBAD c\n" as *u8) } 178 if bf_tan_f64(0x7FF0000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; btg_puts("SPECBAD d\n" as *u8) } 179 if bf_tan_f64(0x4140000000000000) != 0x7FF8000000000000 { sbad = sbad + 1; btg_puts("SPECBAD e\n" as *u8) } // 2^21 refused 180 let below: i64 = bf_tan_f64(0x3FF921FB54442D18) // fl(pi/2) < pi/2: huge + 181 let above: i64 = bf_tan_f64(0x3FF921FB54442D19) // next ulp: crossed pole, - 182 if ((below >> 63) & 1) != 0 { sbad = sbad + 1; btg_puts("SPECBAD f\n" as *u8) } 183 if ((below >> 52) & 0x7FF) < 1073 { sbad = sbad + 1; btg_puts("SPECBAD g\n" as *u8) } // |tan| > 2^50 184 if ((above >> 63) & 1) != 1 { sbad = sbad + 1; btg_puts("SPECBAD h\n" as *u8) } 185 if ((above >> 52) & 0x7FF) < 1073 { sbad = sbad + 1; btg_puts("SPECBAD i\n" as *u8) } 186 let pv: i64 = bf_tan_f64(0x3FE0000000000000) 187 if bf_tan_f64(0xBFE0000000000000) != (pv | (1 << 63)) { sbad = sbad + 1; btg_puts("SPECBAD j\n" as *u8) } 188 btg_puts("BFTAN-SPEC bad=" as *u8); btg_putn(sbad); btg_puts("\n" as *u8) 189 bad = bad + sbad 190 191 if bad == 0 { 192 btg_puts("BFTAN-GATE verdict=GREEN\n" as *u8) 193 return 0 194 } 195 btg_puts("BFTAN-GATE verdict=RED\n" as *u8) 196 if bad > 100 { bad = 100 } 197 return bad 198}