code wiki / (root) / nx_vecmath_gate.nx
`nx_vecmath_gate.nx` source

↩ module page · 344 lines · 18215 B
1// nx_vecmath_gate.nx -- THE GATE for the shared integer linear-algebra foundation (PG2, procgen.plan).
2//
3// GATE LAW APPLIED (CLAUDE.md "EVERY GATE YOU WRITE"):
4//  1. returns gv_verdict so the EXIT CODE carries the verdict -- /api/gate_run derives GREEN/RED from it,
5//     and a gate that returns a bare 0 after printing RED silently blesses every failure it finds.
6//  2. two teeth are named neg-control-* so the L2 census can SEE them.
7//  3. every tooth is a gv_check on the shared counter, so DECLARED == EXECUTED by construction. A hand-rolled
8//     tally can print "passed 22/20"; this cannot.
9//  4. every sweep is bound to its denominator via gv_subjects. A tooth that passes on the empty set is not a
10//     tooth, and these sweeps are the whole representable population, never a sample.
11//  5. every tooth PRINTS ITS MEASURED VALUE. A boolean cannot say why, and both vacuous teeth this estate has
12//     caught were found in the diagnostic dump, never in the verdict vector.
13//  6. the verdict note does NOT recite the teeth. Each gv_check name states its own strength, so adding a
14//     tooth can never make a summary sentence lie.
15//
16// ERROR BARS ARE DERIVED AND THE MEASUREMENT IS PRINTED BESIDE THEM. nx_itrig declares ~0.6 percent amplitude
17// error, which at fx4096 is about 25 units of value and therefore about 25 units of ANGLE once a point is
18// perturbed off the unit circle. The single-term atan form contributes about 20 units. The bars below sit at
19// roughly twice their sum so a tooth fails on a REGRESSION rather than on the known approximation, and the
20// printed maximum makes the bar ratchetable on evidence instead of by taste.
21
22import "nx_syscalls.nx"
23import "nx_gate_verdict.nx"
24import "nx_vecmath.nx"
25
26const VG_I64: i64 = 8
27const VG_V3B: i64 = 24          // 3 * VG_I64
28const VG_V4B: i64 = 32          // 4 * VG_I64
29const VG_M3B: i64 = 72          // 9 * VG_I64
30
31// Derived bars -- see the header. Each is a REGRESSION bar, not a claim of accuracy.
32const VG_ANG_BAR:  i64 = 96     // composite sin-cos-atan2 round-trip, units of rad*4096
33const VG_ACOS_BAR: i64 = 96     // acos(cos(a)) round-trip, same basis
34const VG_UNIT_BAR: i64 = 24     // |normalise(v)| vs VM_ONE: pure fixed-point quantisation, no trig involved
35const VG_ORTH_BAR: i64 = 24     // dot(cross(a,b), a) vs zero
36const VG_M3_BAR:   i64 = 48     // matrix round-trip and rotation-path agreement
37// Polar amplification floor for the acos round-trip: within the cap where sin(theta) < 1/16 the derived
38// bar is capped at 16x the base bar. DERIVED: 1/16 is the first power-of-two floor above the measured
39// worst amplification that still fails on a genuine regression (a dropped quadrant reads ~thousands).
40const VG_POLE_AMP: i64 = 16
41
42// Sweep extents. Angles sweep the FULL representable circle at unit resolution: this is the population,
43// not a sample of it.
44const VG_ISQRT_N:  i64 = 3000
45const VG_FLOOR_N:  i64 = 20000
46const VG_ACOS_N:   i64 = 4096   // every representable cosine in [0, VM_ONE]
47
48func vg_dev(a: i64, b: i64) -> i64 { if a > b { return a - b } return b - a }
49
50// Angular difference on the circle, so a wrap from +PI to -PI reads as small, not as 2*PI.
51func vg_angdev(a: i64, b: i64) -> i64 {
  var d: i64 = a - b
  while d >  vm_pi() { d = d - vm_tau() }
  while d < 0 - vm_pi() { d = d + vm_tau() }
  return vm_abs(d)
56}
57
58func main() -> i64 {
  let c: *i64 = gv_ctr()
  gv_head("NX-VECMATH-GATE -- shared integer linear algebra (PG2)" as *u8)
61
  let one: i64 = vm_one()
63
  // ---- T1: the duplicated angle constants must agree with the incumbent they restate --------------
  // VM_PI and VM_PI2 are a DELIBERATE duplicate of nx_itrig's constants because const scope does not
  // cross the import. This tooth is what makes the duplicate guarded rather than silent.
  let sin_pi:  i64 = vm_sin(vm_pi())
  let sin_pi2: i64 = vm_sin(vm_pi2())
  gv_puts("  pi_check sin(PI)=" as *u8); gv_num(sin_pi)
  gv_puts(" sin(PI/2)=" as *u8); gv_num(sin_pi2)
  gv_puts(" expect ~0 and ~" as *u8); gv_num(one); gv_puts("\n" as *u8)
  gv_check("pi-agrees-with-itrig" as *u8,
           vm_abs(sin_pi) < VG_ANG_BAR && vg_dev(sin_pi2, one) < VG_ANG_BAR, c)
74
  // ---- T2/T3: isqrt, full sweep ------------------------------------------------------------------
  gv_subjects("isqrt-perfect-squares" as *u8, VG_ISQRT_N, c)
  var bad_sq: i64 = 0
  var n: i64 = 0
  while n < VG_ISQRT_N {
      if vm_isqrt(n * n) != n { bad_sq = bad_sq + 1 }
      n = n + 1
  }
  gv_puts("  isqrt exact_failures=" as *u8); gv_num(bad_sq)
  gv_puts(" of " as *u8); gv_num(VG_ISQRT_N); gv_puts("\n" as *u8)
  gv_check("isqrt-exact-on-every-perfect-square" as *u8, bad_sq == 0, c)
86
  // The DEFINING property, stronger than agreeing with any incumbent: r*r <= v < (r+1)*(r+1).
  gv_subjects("isqrt-floor-property" as *u8, VG_FLOOR_N, c)
  var bad_fl: i64 = 0
  var v: i64 = 1
  while v <= VG_FLOOR_N {
      let r: i64 = vm_isqrt(v)
      if r * r > v { bad_fl = bad_fl + 1 }
      if (r + 1) * (r + 1) <= v { bad_fl = bad_fl + 1 }
      v = v + 1
  }
  gv_puts("  isqrt floor_violations=" as *u8); gv_num(bad_fl)
  gv_puts(" of " as *u8); gv_num(VG_FLOOR_N); gv_puts("\n" as *u8)
  gv_check("isqrt-is-exact-floor-not-approximate" as *u8, bad_fl == 0, c)
100
  // ---- T4: atan2 round-trip over the FULL circle -------------------------------------------------
  // Measures the COMPOSITE sin-cos-atan2 error, which is what every consumer actually experiences.
  let sweep: i64 = vm_tau()
  gv_subjects("atan2-roundtrip-full-circle" as *u8, sweep, c)
  var worst_ang: i64 = 0
  var worst_at: i64 = 0
  var a: i64 = 0 - vm_pi() + 1
  var seen: i64 = 0
  while a <= vm_pi() {
      let s: i64 = vm_sin(a)
      let co: i64 = vm_cos(a)
      let back: i64 = vm_atan2(s, co)
      let d: i64 = vg_angdev(back, a)
      if d > worst_ang { worst_ang = d; worst_at = a }
      seen = seen + 1
      a = a + 1
  }
  gv_puts("  atan2 worst_err=" as *u8); gv_num(worst_ang)
  gv_puts(" at_angle=" as *u8); gv_num(worst_at)
  gv_puts(" bar=" as *u8); gv_num(VG_ANG_BAR)
  gv_puts(" sampled=" as *u8); gv_num(seen); gv_puts("\n" as *u8)
  // Bound to its denominator: a sweep that examined nothing must not pass.
  gv_check("atan2-roundtrip-within-derived-bar-over-full-circle" as *u8,
           seen > 0 && worst_ang < VG_ANG_BAR, c)
125
  // ---- T5: quadrant signs. A single-quadrant atan is a common silent defect. ----------------------
  let q1: i64 = vm_atan2(one, one)
  let q2: i64 = vm_atan2(one, 0 - one)
  let q3: i64 = vm_atan2(0 - one, 0 - one)
  let q4: i64 = vm_atan2(0 - one, one)
  gv_puts("  quadrants q1=" as *u8); gv_num(q1); gv_puts(" q2=" as *u8); gv_num(q2)
  gv_puts(" q3=" as *u8); gv_num(q3); gv_puts(" q4=" as *u8); gv_num(q4); gv_puts("\n" as *u8)
  gv_check("atan2-covers-all-four-quadrants" as *u8,
           q1 > 0 && q2 > 0 && q3 < 0 && q4 < 0 && q2 > q1 && q3 < q4, c)
135
  // ---- T6: acos, exact algebra first, then the round-trip with a bar DERIVED per sample ------------
  // FIRST RUN 2026-08-24 read worst_err=334 against a constant bar of 96 and went RED -- and the bar was
  // the defect, not the library. d(acos)/dc = 1/sin(theta): near the poles a 25-unit error in vm_cos's
  // Taylor output is amplified up to 70x into angle. A constant bar assumed unit slope. So the round-trip
  // bar is now the amplification itself, derived per sample, floored at 16x within the polar cap where
  // sin(theta) < 1/16 -- and the exact identities below carry no incumbent error at all.
  let acos_top: i64 = vm_acos(one)
  let acos_mid: i64 = vm_acos(0)
  let acos_bot: i64 = vm_acos(0 - one)
  gv_puts("  acos(1)=" as *u8); gv_num(acos_top); gv_puts(" acos(0)=" as *u8); gv_num(acos_mid)
  gv_puts(" acos(-1)=" as *u8); gv_num(acos_bot)
  gv_puts(" expect 0 " as *u8); gv_num(vm_pi2()); gv_puts(" " as *u8); gv_num(vm_pi()); gv_puts("\n" as *u8)
  gv_check("acos-hits-its-three-exact-anchors" as *u8,
           acos_top == 0 && vg_dev(acos_mid, vm_pi2()) <= 1 && vg_dev(acos_bot, vm_pi()) <= 1, c)
150
  // Monotone decreasing over EVERY representable cosine -- the full population, exact, no trig error.
  let acos_pop: i64 = 2 * one + 1
  gv_subjects("acos-monotone-full-domain" as *u8, acos_pop, c)
  var mono_viol: i64 = 0
  var prev_ac: i64 = vm_acos(0 - one)
  var cc: i64 = 0 - one + 1
  while cc <= one {
      let cur: i64 = vm_acos(cc)
      if cur > prev_ac { mono_viol = mono_viol + 1 }
      prev_ac = cur
      cc = cc + 1
  }
  gv_puts("  acos monotone_violations=" as *u8); gv_num(mono_viol); gv_puts(" of " as *u8); gv_num(acos_pop); gv_puts("\n" as *u8)
  gv_check("acos-is-monotone-decreasing-over-every-representable-input" as *u8, mono_viol == 0, c)
165
  gv_subjects("acos-roundtrip" as *u8, VG_ACOS_N, c)
  var worst_ratio: i64 = 0        // err * sin(theta) / bar, permil: 1000 = exactly at the derived bar
  var worst_ac: i64 = 0
  var worst_ac_at: i64 = 0
  var ac_seen: i64 = 0
  var ac_viol: i64 = 0
  var k: i64 = 0
  while k < VG_ACOS_N {
      let ang: i64 = vm_pi() * k / VG_ACOS_N
      let cv: i64 = vm_cos(ang)
      let back2: i64 = vm_acos(cv)
      let d2: i64 = vg_dev(back2, ang)
      var sn: i64 = vm_sin(ang)
      if sn < one / VG_POLE_AMP { sn = one / VG_POLE_AMP }
      let allowed: i64 = VG_ACOS_BAR * one / sn
      if d2 >= allowed { ac_viol = ac_viol + 1 }
      let ratio: i64 = d2 * 1000 / allowed
      if ratio > worst_ratio { worst_ratio = ratio }
      if d2 > worst_ac { worst_ac = d2; worst_ac_at = ang }
      ac_seen = ac_seen + 1
      k = k + 1
  }
  gv_puts("  acos worst_err=" as *u8); gv_num(worst_ac); gv_puts(" at_angle=" as *u8); gv_num(worst_ac_at)
  gv_puts(" worst_err_over_derived_bar_permil=" as *u8); gv_num(worst_ratio)
  gv_puts(" violations=" as *u8); gv_num(ac_viol)
  gv_puts(" base_bar=" as *u8); gv_num(VG_ACOS_BAR); gv_puts(" pole_amp_floor=" as *u8); gv_num(VG_POLE_AMP)
  gv_puts(" sampled=" as *u8); gv_num(ac_seen); gv_puts("\n" as *u8)
  gv_check("acos-roundtrip-within-the-per-sample-amplification-derived-bar" as *u8,
           ac_seen > 0 && ac_viol == 0, c)
195
  // ---- vec3 ---------------------------------------------------------------------------------------
  let va: *i64 = sys_mmap(VG_V3B) as *i64
  let vb: *i64 = sys_mmap(VG_V3B) as *i64
  let vc: *i64 = sys_mmap(VG_V3B) as *i64
  let vd: *i64 = sys_mmap(VG_V3B) as *i64
201
  gv_subjects("v3-cross-orthogonality" as *u8, VG_ACOS_N, c)
  var worst_orth: i64 = 0
  var worst_unit: i64 = 0
  var vseen: i64 = 0
  var t: i64 = 0
  while t < VG_ACOS_N {
      let ang2: i64 = vm_tau() * t / VG_ACOS_N
      vm_v3_set(va, vm_cos(ang2), vm_sin(ang2), one / 2)
      vm_v3_set(vb, 0 - vm_sin(ang2), vm_cos(ang2), 0 - one / 3)
      vm_v3_cross(vc, va, vb)
      let da: i64 = vm_abs(vm_v3_dot(vc, va))
      let db: i64 = vm_abs(vm_v3_dot(vc, vb))
      if da > worst_orth { worst_orth = da }
      if db > worst_orth { worst_orth = db }
      if vm_v3_norm(vd, va) == 1 {
          let ul: i64 = vg_dev(vm_v3_len(vd), one)
          if ul > worst_unit { worst_unit = ul }
      }
      vseen = vseen + 1
      t = t + 1
  }
  gv_puts("  v3 worst_orth=" as *u8); gv_num(worst_orth)
  gv_puts(" worst_unit_dev=" as *u8); gv_num(worst_unit)
  gv_puts(" sampled=" as *u8); gv_num(vseen); gv_puts("\n" as *u8)
  gv_check("v3-cross-is-orthogonal-to-both-operands" as *u8,
           vseen > 0 && worst_orth < VG_ORTH_BAR, c)
  gv_check("v3-normalise-yields-unit-length" as *u8,
           vseen > 0 && worst_unit < VG_UNIT_BAR, c)
230
  // neg-control: the zero vector has no direction. Returning 1 here, or writing garbage, is the defect.
  vm_v3_set(va, 0, 0, 0)
  vm_v3_set(vd, 12345, 12345, 12345)
  let zr: i64 = vm_v3_norm(vd, va)
  gv_puts("  neg-control zero_norm_ret=" as *u8); gv_num(zr)
  gv_puts(" out=" as *u8); gv_num(vd[0]); gv_puts("," as *u8); gv_num(vd[1])
  gv_puts("," as *u8); gv_num(vd[2]); gv_puts("\n" as *u8)
  gv_check("neg-control-v3-normalise-refuses-zero-vector-and-zeroes-output" as *u8,
           zr == 0 && vd[0] == 0 && vd[1] == 0 && vd[2] == 0, c)
240
  // ---- quaternion ---------------------------------------------------------------------------------
  let qa: *i64 = sys_mmap(VG_V4B) as *i64
  let qb: *i64 = sys_mmap(VG_V4B) as *i64
  let qr: *i64 = sys_mmap(VG_V4B) as *i64
  let qi: *i64 = sys_mmap(VG_V4B) as *i64
246
  vm_v3_set(va, 0, one, 0)
  vm_q_from_axis(qa, va, vm_pi2())          // 90 deg about +Y
  vm_q_ident(qi)
250
  vm_q_mul(qr, qa, qi)
  let dmul: i64 = vg_dev(qr[0], qa[0]) + vg_dev(qr[1], qa[1]) + vg_dev(qr[2], qa[2]) + vg_dev(qr[3], qa[3])
  gv_puts("  q*ident dev=" as *u8); gv_num(dmul); gv_puts("\n" as *u8)
  gv_check("q-mul-identity-is-identity" as *u8, dmul < VG_UNIT_BAR, c)
255
  vm_q_conj(qb, qa)
  vm_q_mul(qr, qa, qb)
  let dconj: i64 = vg_dev(qr[3], one) + vm_abs(qr[0]) + vm_abs(qr[1]) + vm_abs(qr[2])
  gv_puts("  q*conj dev_from_ident=" as *u8); gv_num(dconj); gv_puts("\n" as *u8)
  gv_check("q-times-its-conjugate-is-identity" as *u8, dconj < VG_M3_BAR, c)
261
  // The two rotation paths must agree. If they diverge, one of them is wrong and neither is trustworthy.
  let m1: *i64 = sys_mmap(VG_M3B) as *i64
  let m2: *i64 = sys_mmap(VG_M3B) as *i64
  vm_v3_set(vb, one, 0, 0)
  vm_q_rotate_v3(vc, qa, vb)
  vm_q_to_m3(m1, qa)
  vm_m3_mulv(vd, m1, vb)
  let dpath: i64 = vg_dev(vc[0], vd[0]) + vg_dev(vc[1], vd[1]) + vg_dev(vc[2], vd[2])
  gv_puts("  rotate_v3=" as *u8); gv_num(vc[0]); gv_puts("," as *u8); gv_num(vc[1])
  gv_puts("," as *u8); gv_num(vc[2])
  gv_puts("  m3*v=" as *u8); gv_num(vd[0]); gv_puts("," as *u8); gv_num(vd[1])
  gv_puts("," as *u8); gv_num(vd[2])
  gv_puts("  dev=" as *u8); gv_num(dpath); gv_puts("\n" as *u8)
  gv_check("quaternion-and-matrix-rotation-paths-agree" as *u8, dpath < VG_M3_BAR, c)
276
  // ---- slerp --------------------------------------------------------------------------------------
  vm_v3_set(va, 0, one, 0)
  vm_q_from_axis(qb, va, vm_pi2())
  vm_q_slerp(qr, qi, qb, 0)
  let d_lo: i64 = vg_dev(qr[3], one) + vm_abs(qr[0]) + vm_abs(qr[1]) + vm_abs(qr[2])
  vm_q_slerp(qr, qi, qb, one)
  let d_hi: i64 = vg_dev(qr[0], qb[0]) + vg_dev(qr[1], qb[1]) + vg_dev(qr[2], qb[2]) + vg_dev(qr[3], qb[3])
  gv_puts("  slerp t=0 dev=" as *u8); gv_num(d_lo)
  gv_puts(" t=ONE dev=" as *u8); gv_num(d_hi); gv_puts("\n" as *u8)
  gv_check("slerp-endpoints-return-the-endpoints" as *u8,
           d_lo < VG_M3_BAR && d_hi < VG_M3_BAR, c)
288
  // ANTI-VACUITY. slerp falls back to normalised lerp for near-parallel inputs. If that fallback ran for
  // EVERY input the spherical arithmetic would be dead code and every tooth above would still pass. For a
  // 90-degree pair the spherical midpoint and the lerp midpoint are genuinely different, so this tooth
  // fails if the fallback has swallowed the real path. A trivial lerp-only implementation CANNOT pass it.
  vm_q_slerp(qr, qi, qb, one / 2)
  let lx: i64 = (qi[0] + qb[0]) / 2
  let ly: i64 = (qi[1] + qb[1]) / 2
  let lz: i64 = (qi[2] + qb[2]) / 2
  let lw: i64 = (qi[3] + qb[3]) / 2
  let dlerp: i64 = vg_dev(qr[0], lx) + vg_dev(qr[1], ly) + vg_dev(qr[2], lz) + vg_dev(qr[3], lw)
  let mid_len: i64 = vm_q_len(qr)
  gv_puts("  slerp_mid vs raw_lerp_mid dev=" as *u8); gv_num(dlerp)
  gv_puts(" mid_len=" as *u8); gv_num(mid_len)
  gv_puts(" expect dev>0 and len~" as *u8); gv_num(one); gv_puts("\n" as *u8)
  gv_check("anti-vacuity-slerp-spherical-branch-actually-runs" as *u8,
           dlerp > 0 && vg_dev(mid_len, one) < VG_UNIT_BAR, c)
305
  // neg-control: an axis with no direction cannot define a rotation.
  vm_v3_set(va, 0, 0, 0)
  let far: i64 = vm_q_from_axis(qr, va, vm_pi2())
  gv_puts("  neg-control zero_axis_ret=" as *u8); gv_num(far)
  gv_puts(" q=" as *u8); gv_num(qr[0]); gv_puts("," as *u8); gv_num(qr[1])
  gv_puts("," as *u8); gv_num(qr[2]); gv_puts("," as *u8); gv_num(qr[3]); gv_puts("\n" as *u8)
  gv_check("neg-control-q-from-zero-axis-refuses-and-returns-identity" as *u8,
           far == 0 && qr[3] == one && qr[0] == 0 && qr[1] == 0 && qr[2] == 0, c)
314
  // ---- mat3 ---------------------------------------------------------------------------------------
  vm_m3_ident(m2)
  vm_q_to_m3(m1, qa)
  let m1b: i64 = m1[0]; let m1c: i64 = m1[4]; let m1d: i64 = m1[8]
  vm_m3_mul(m1, m1, m2)
  let dident: i64 = vg_dev(m1[0], m1b) + vg_dev(m1[4], m1c) + vg_dev(m1[8], m1d)
  gv_puts("  m3*I dev=" as *u8); gv_num(dident); gv_puts("\n" as *u8)
  gv_check("m3-multiply-by-identity-is-a-no-op-and-alias-safe" as *u8, dident < VG_M3_BAR, c)
323
  vm_q_to_m3(m1, qa)
  let orig0: i64 = m1[1]; let orig1: i64 = m1[5]
  vm_m3_transpose(m2, m1)
  vm_m3_transpose(m2, m2)
  let dtr: i64 = vg_dev(m2[1], orig0) + vg_dev(m2[5], orig1)
  gv_puts("  m3 transpose-twice dev=" as *u8); gv_num(dtr); gv_puts("\n" as *u8)
  gv_check("m3-transpose-is-an-involution-and-alias-safe" as *u8, dtr == 0, c)
331
  // The covariance builder: the call nx_gsplat could not make before this lib existed.
  vm_q_scale_to_m3(m1, qi, one * 2, one * 3, one * 4)
  gv_puts("  q_scale_to_m3 diag=" as *u8); gv_num(m1[0]); gv_puts("," as *u8); gv_num(m1[4])
  gv_puts("," as *u8); gv_num(m1[8]); gv_puts(" expect " as *u8); gv_num(one*2)
  gv_puts("," as *u8); gv_num(one*3); gv_puts("," as *u8); gv_num(one*4); gv_puts("\n" as *u8)
  gv_check("q-scale-to-m3-applies-all-three-axis-scales" as *u8,
           vg_dev(m1[0], one*2) < VG_M3_BAR &&
           vg_dev(m1[4], one*3) < VG_M3_BAR &&
           vg_dev(m1[8], one*4) < VG_M3_BAR, c)
341
  return gv_verdict("nx_vecmath_gate" as *u8, c,
      "Shared integer linear algebra. Bars are REGRESSION bars derived from the incumbent's declared error, not accuracy claims; each tooth prints its measured value so the bars stay ratchetable on evidence." as *u8)
344}
nx_vecmath_gate.nx source

`nx_vecmath_gate.nx` source
