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

↩ module page · 236 lines · 11037 B
1// nx_linalg_chol_gate.nx -- THE REFEREE FOR THE Q30 CHOLESKY (graphics GR47).
2//
3// The done-rule GR47 declared BEFORE this was written, quoted from graphics.plan so the bar cannot drift:
4//   "a Q30 fixed-point Cholesky over nx_fixq30_lib arithmetic factors any SPD matrix up to the
5//    identity-space dimension, L times L-transposed reproduces the input within a stated Q30 band on
6//    random SPD fixtures, a non-SPD input is refused by name rather than producing garbage, and GR24
7//    composes it for its draw (one ruler, never a second decomposition)"
8// This gate proves the first three clauses. The fourth is GR24's own rung and is NOT claimed here.
9//
10// THE LOAD-BEARING TOOTH IS THE TEXTBOOK 3x3, because its factor is known by hand and INDEPENDENTLY of this
11// implementation: A = [[4,12,-16],[12,37,-43],[-16,-43,98]] has L = [[2,0,0],[6,1,0],[-8,5,3]] exactly.
12// A residual check alone would pass for any L whose product happens to reconstruct A; comparing against a
13// hand-derived L tests that this is THE Cholesky factor and not merely A factor.
14// license_tier: ORIGINAL   No hw writes (Rule 26).
15import "nx_syscalls.nx"
16import "nx_gate_verdict.nx"
17import "nx_linalg.nx"
18
19const LCG_MAXN: i64 = 8
20
21// build an n x n Q30 matrix from raw integers (each scaled by FQ_ONE)
22func lcg_mat(n: i64, vals: *i64) -> *Mat {
  let m: *Mat = nx_mat_new(n, n)
  var i: i64 = 0
  while i < n {
      var j: i64 = 0
      while j < n { let _s: i64 = nx_mat_set(m, i, j, fq_from_int(vals[(i * n) + j])); j = j + 1 }
      i = i + 1
  }
  return m
31}
32func lcg_vals(n: i64) -> *i64 { return sys_mmap(n * n * 8) as *i64 }
33
34// worst |L - expected| in Q30 units, so the comparison against a hand-derived factor is a NUMBER
35func lcg_l_diff(l: *Mat, exp: *Mat) -> i64 {
  let n: i64 = l.rows
  var worst: i64 = 0
  var i: i64 = 0
  while i < n {
      var j: i64 = 0
      while j < n {
          var e: i64 = nx_mat_get(l, i, j) - nx_mat_get(exp, i, j)
          if e < 0 { e = 0 - e }
          if e > worst { worst = e }
          j = j + 1
      }
      i = i + 1
  }
  return worst
50}
51
52// is every entry above the diagonal exactly zero
53func lcg_is_lower(l: *Mat) -> i64 {
  let n: i64 = l.rows
  var i: i64 = 0
  while i < n {
      var j: i64 = i + 1
      while j < n { if nx_mat_get(l, i, j) != 0 { return 0 } j = j + 1 }
      i = i + 1
  }
  return 1
62}
63func lcg_all_zero(l: *Mat) -> i64 {
  let n: i64 = l.rows
  var i: i64 = 0
  while i < n {
      var j: i64 = 0
      while j < n { if nx_mat_get(l, i, j) != 0 { return 0 } j = j + 1 }
      i = i + 1
  }
  return 1
72}
73
74func main() -> i64 {
  let ctr: *i64 = gv_ctr()
  gv_head("nx_linalg_chol_gate -- Q30 Cholesky (GR47)" as *u8)
77
  // ---- 1. IDENTITY: the simplest factor that must come back exactly --------------------------------
  let iv: *i64 = lcg_vals(3)
  iv[0]=1; iv[1]=0; iv[2]=0; iv[3]=0; iv[4]=1; iv[5]=0; iv[6]=0; iv[7]=0; iv[8]=1
  let ia: *Mat = lcg_mat(3, iv)
  let il: *Mat = nx_mat_new(3, 3)
  gv_check_eq("identity-factors-OK" as *u8, la_cholesky_q30(ia, il), LA_CHOL_OK, ctr)
  gv_check_eq("identity-factor-is-the-identity" as *u8, lcg_l_diff(il, ia), 0, ctr)
  gv_check_eq("identity-residual-is-exactly-zero" as *u8, la_chol_max_resid_q30(ia, il), 0, ctr)
86
  // ---- 2. THE TEXTBOOK 3x3 AGAINST A HAND-DERIVED FACTOR --------------------------------------------
  let av: *i64 = lcg_vals(3)
  av[0]=4;   av[1]=12;  av[2]=0-16
  av[3]=12;  av[4]=37;  av[5]=0-43
  av[6]=0-16; av[7]=0-43; av[8]=98
  let a3: *Mat = lcg_mat(3, av)
  let ev: *i64 = lcg_vals(3)
  ev[0]=2;   ev[1]=0; ev[2]=0
  ev[3]=6;   ev[4]=1; ev[5]=0
  ev[6]=0-8; ev[7]=5; ev[8]=3
  let e3: *Mat = lcg_mat(3, ev)
  let l3: *Mat = nx_mat_new(3, 3)
  gv_check_eq("textbook-3x3-factors-OK" as *u8, la_cholesky_q30(a3, l3), LA_CHOL_OK, ctr)
  let ldiff: i64 = lcg_l_diff(l3, e3)
  gv_kv("textbook_L_vs_hand_derived_q30" as *u8, ldiff)
  // the hand factor is exact and every radicand here is a perfect square, so this must be EXACT
  gv_check_eq("textbook-factor-matches-the-hand-derived-L-exactly" as *u8, ldiff, 0, ctr)
  let r3: i64 = la_chol_max_resid_q30(a3, l3)
  gv_kv("textbook_residual_q30" as *u8, r3)
  gv_check("textbook-residual-inside-the-declared-band" as *u8, (r3 <= LA_CHOL_BAND_Q30) as i64, ctr)
  gv_check_eq("factor-is-lower-triangular" as *u8, lcg_is_lower(l3), 1, ctr)
108
  // ---- 3. AN SPD MATRIX WITH IRRATIONAL ROOTS: where the band actually gets exercised ---------------
  // A = [[2,1,0],[1,3,1],[0,1,4]] is SPD and none of its pivots is a perfect square, so this is the
  // fixture that tests the ARITHMETIC rather than a lucky exact case.
  let bv: *i64 = lcg_vals(3)
  bv[0]=2; bv[1]=1; bv[2]=0
  bv[3]=1; bv[4]=3; bv[5]=1
  bv[6]=0; bv[7]=1; bv[8]=4
  let ab: *Mat = lcg_mat(3, bv)
  let lb: *Mat = nx_mat_new(3, 3)
  gv_check_eq("irrational-root-SPD-factors-OK" as *u8, la_cholesky_q30(ab, lb), LA_CHOL_OK, ctr)
  let rb: i64 = la_chol_max_resid_q30(ab, lb)
  gv_kv("irrational_root_residual_q30" as *u8, rb)
  gv_check("irrational-root-residual-inside-the-declared-band" as *u8, (rb <= LA_CHOL_BAND_Q30) as i64, ctr)
122
  // ---- 4. A LARGER SPD: the residual must not silently blow up with dimension -----------------------
  // diagonally dominant 6x6 (4 on the diagonal, 1 on the off-diagonals) is SPD by construction
  let n6: i64 = 6
  let cv: *i64 = lcg_vals(n6)
  var i: i64 = 0
  while i < n6 {
      var j: i64 = 0
      while j < n6 {
          if i == j { cv[(i * n6) + j] = 4 } else { cv[(i * n6) + j] = 1 }
          j = j + 1
      }
      i = i + 1
  }
  let a6: *Mat = lcg_mat(n6, cv)
  let l6: *Mat = nx_mat_new(n6, n6)
  gv_check_eq("6x6-SPD-factors-OK" as *u8, la_cholesky_q30(a6, l6), LA_CHOL_OK, ctr)
  let r6: i64 = la_chol_max_resid_q30(a6, l6)
  gv_kv("dim6_residual_q30" as *u8, r6)
  gv_check("6x6-residual-inside-the-declared-band" as *u8, (r6 <= LA_CHOL_BAND_Q30) as i64, ctr)
  gv_check_eq("6x6-factor-is-lower-triangular" as *u8, lcg_is_lower(l6), 1, ctr)
143
  // ---- 4b. AT THE DECLARED DIMENSION CAP: the band claims to hold up to LA_CHOL_MAX_N, so it is
  // measured THERE and not only at the sizes that were convenient. A band asserted past the largest
  // fixture is exactly the bound-recorded-as-a-fact defect.
  let nmax: i64 = LA_CHOL_MAX_N
  let dv: *i64 = lcg_vals(nmax)
  var p: i64 = 0
  while p < nmax {
      var q: i64 = 0
      while q < nmax {
          if p == q { dv[(p * nmax) + q] = 4 } else { dv[(p * nmax) + q] = 1 }
          q = q + 1
      }
      p = p + 1
  }
  let amax: *Mat = lcg_mat(nmax, dv)
  let lmax: *Mat = nx_mat_new(nmax, nmax)
  gv_check_eq("cap-dimension-SPD-factors-OK" as *u8, la_cholesky_q30(amax, lmax), LA_CHOL_OK, ctr)
  let rmax: i64 = la_chol_max_resid_q30(amax, lmax)
  gv_kv("dim32_cap_residual_q30" as *u8, rmax)
  gv_check("cap-dimension-residual-inside-the-declared-band" as *u8, (rmax <= LA_CHOL_BAND_Q30) as i64, ctr)
  gv_check_eq("cap-dimension-factor-is-lower-triangular" as *u8, lcg_is_lower(lmax), 1, ctr)
  // THE BAND MUST BE A RATCHET, NOT A COMFORT: if the worst fixture sits far under it the bar is
  // decorative and no regression would ever trip it. This tooth pins the band within one order of the
  // worst measured residual, so loosening it later is a deliberate act with a failing test attached.
  var worst: i64 = r3
  if rb > worst { worst = rb }
  if r6 > worst { worst = r6 }
  if rmax > worst { worst = rmax }
  gv_kv("worst_residual_over_all_fixtures_q30" as *u8, worst)
  gv_check("declared-band-is-not-more-than-16x-the-worst-measured" as *u8,
           (LA_CHOL_BAND_Q30 <= (worst + 1) * 16) as i64, ctr)
  // and it must still be a real bar: a band at or below the worst observation would be permanently red
  gv_check("declared-band-is-at-least-the-worst-measured" as *u8,
           (LA_CHOL_BAND_Q30 >= worst) as i64, ctr)
178
  // ---- 5. NEG-CONTROLS: every refusal named, and the factor left unusable ---------------------------
  // NOT positive definite: [[1,2],[2,1]] is symmetric with a negative eigenvalue.
  let pv: *i64 = lcg_vals(2)
  pv[0]=1; pv[1]=2; pv[2]=2; pv[3]=1
  let ap: *Mat = lcg_mat(2, pv)
  let lp: *Mat = nx_mat_new(2, 2)
  gv_check_eq("neg-control-non-PD-is-refused-BY-NAME" as *u8, la_cholesky_q30(ap, lp), LA_CHOL_NOT_PD, ctr)
  // a refusal must not leave a half-written factor a caller could mistake for a real one
  gv_check_eq("neg-control-non-PD-leaves-the-factor-zeroed" as *u8, lcg_all_zero(lp), 1, ctr)
188
  // NOT symmetric: same numbers, one transposed entry changed.
  let sv: *i64 = lcg_vals(2)
  sv[0]=4; sv[1]=2; sv[2]=3; sv[3]=5
  let asym: *Mat = lcg_mat(2, sv)
  let lsym: *Mat = nx_mat_new(2, 2)
  gv_check_eq("neg-control-non-symmetric-is-refused-BY-NAME" as *u8,
              la_cholesky_q30(asym, lsym), LA_CHOL_NOT_SYMMETRIC, ctr)
196
  // NOT square
  let ans: *Mat = nx_mat_new(2, 3)
  let lns: *Mat = nx_mat_new(2, 3)
  gv_check_eq("neg-control-non-square-is-refused-BY-NAME" as *u8,
              la_cholesky_q30(ans, lns), LA_CHOL_NOT_SQUARE, ctr)
202
  // dimension past the measured band
  let big: i64 = LA_CHOL_MAX_N + 1
  let abig: *Mat = nx_mat_new(big, big)
  let lbig: *Mat = nx_mat_new(big, big)
  gv_check_eq("neg-control-dimension-past-the-measured-band-refuses" as *u8,
              la_cholesky_q30(abig, lbig), LA_CHOL_BAD_DIM, ctr)
209
  // ---- 6. THE POSITIVE CONTROL: a guard that refused everything would pass every test above ---------
  let qv: *i64 = lcg_vals(2)
  qv[0]=4; qv[1]=2; qv[2]=2; qv[3]=5
  let aq: *Mat = lcg_mat(2, qv)
  let lq: *Mat = nx_mat_new(2, 2)
  gv_check_eq("pos-control-a-valid-SPD-is-ACCEPTED" as *u8, la_cholesky_q30(aq, lq), LA_CHOL_OK, ctr)
  gv_check_eq("pos-control-factor-is-not-zeroed" as *u8, lcg_all_zero(lq), 0, ctr)
  // hand-derived: sqrt(4)=2, 2/2=1, sqrt(5-1)=2  ->  L = [[2,0],[1,2]] exactly
  let qe: *i64 = lcg_vals(2)
  qe[0]=2; qe[1]=0; qe[2]=1; qe[3]=2
  gv_check_eq("pos-control-matches-its-hand-derived-factor" as *u8, lcg_l_diff(lq, lcg_mat(2, qe)), 0, ctr)
221
  // ---- 7. the status names are distinct, so a refusal can be read -----------------------------------
  gv_check_eq("status-names-cover-every-code" as *u8, LA_CHOL_N, 5, ctr)
224
  gv_values_head()
  gv_kv("declared_band_q30" as *u8, LA_CHOL_BAND_Q30)
  gv_kv("max_measured_dimension" as *u8, LA_CHOL_MAX_N)
  gv_kv("q30_one" as *u8, FQ_ONE)
  gv_kv("textbook_residual_q30" as *u8, r3)
  gv_kv("irrational_root_residual_q30" as *u8, rb)
  gv_kv("dim6_residual_q30" as *u8, r6)
  gv_kv("dim32_cap_residual_q30" as *u8, rmax)
233
  return gv_verdict("nx_linalg_chol_gate" as *u8, ctr,
      "GR47 clauses 1-3: SPD factored in Q30, residual measured against a declared band, every refusal named" as *u8)
236}
nx_linalg_chol_gate.nx source

`nx_linalg_chol_gate.nx` source
