analysis.mac

/* =====================================================================
   analysis.mac

   Maxima analysis and proofs for:
   "Emergent Fractional Dimensionality and Spinor-Like Holonomy
    in Multi-Sheeted n-gon Tilings"

   This file generalizes the original pentagon-specific analysis to an
   arbitrary regular n-gon (n >= 3).  The pentagon (n = 5) case is
   preserved as a verification benchmark and is run automatically at
   the end of the file.

   All computations are performed in exact arithmetic when possible.
   For the n = 5 case, the substrate is the real quadratic field
   Q(sqrt(5)) (equivalently the maximal real subfield of Q[zeta_5]).
   For general n, the substrate is the maximal real subfield of the
   cyclotomic field Q[zeta_n], i.e. Q(zeta_n + zeta_n^{-1}) of
   degree phi(n)/2 over Q (for n >= 3).

   Sections:
     0.  PARAMETERIZATION  (set N_GON here; everything else follows)
     1.  Golden ratio phi and the field Q(sqrt(5))      [pentagon-specific]
     2.  Regular n-gon: angles, deficits, excesses      [generalized]
     3.  n-th roots of unity and the cyclotomic field   [generalized]
     4.  Exact rotations in 2D and basis vectors        [generalized]
     5.  Holonomy, sheet transitions, spinor cover      [generalized]
     6.  Emergent dimensions and Alexander-Orbach       [n-independent]
     7.  Cut-and-project: ambient Z^n lattice           [generalized]
     8.  Cellular automaton: regularity and rule count  [generalized]
     9.  n-gon geometry: diagonals, areas, vertices     [generalized]
    10.  Q(sqrt(5))-arithmetic and units                [pentagon-specific]
    11.  Dihedral group D_n and n-gon symmetries        [generalized]
    12.  Volume growth and dimensional flow             [n-independent]
    13.  Graph Laplacian: star and cycle C_n            [generalized]
    14.  Self-tests and summary
   ===================================================================== */

/* Make all output verbose-friendly. */
if not(?boundp('display2d)) then display2d : false$
/* Load string-handling package (provides pad_right, sconcat utilities). */
load("stringproc")$
/* Load additional packages for enhanced analysis. */
load("eigen")$

print("============================================================")$
print("analysis.mac : Multi-Sheeted n-gon Tiling Analysis")$
print("============================================================")$

/* ---------------------------------------------------------------------
    SECTION 0.  RANGE PARAMETERIZATION
   ---------------------------------------------------------------------

    Three parameters below define the sweep range.  Every n-gon-dependent
    section (2-9, 11-13) is executed inside a loop over n = N_GON_MIN,
    N_GON_MIN + N_GON_STEP, ..., N_GON_MAX.  Pentagon-specific symbolic
    checks (Sections 1, 10) are still executed unconditionally once as a
    benchmark before the loop.

   Quick orientation:

     n = 3  : equilateral triangle (6 around vertex tile flat)
     n = 4  : square              (4 around vertex tile flat)
     n = 5  : regular pentagon    (3 leave deficit pi/5; original paper)
     n = 6  : regular hexagon     (3 tile flat; honeycomb)
     n = 7  : regular heptagon    (3 yield negative-curvature surplus)
     n = 8  : regular octagon
     n = 12 : regular dodecagon   (rich cyclotomic structure)

    Floor( 2*pi / interior_angle ) tells how many n-gons fit at a vertex
    in flat geometry, and the residual angle gives the deficit/excess
    used in the holonomy section.

    RANGE CONTROLS  (edit these three lines):
   --------------------------------------------------------------------- */

N_GON_MIN  : 3  $   /* first n-gon in sweep  (integer >= 3)       */
N_GON_MAX  : 12 $   /* last  n-gon in sweep  (integer >= N_GON_MIN) */
N_GON_STEP : 1  $   /* step size             (positive integer)    */
/* ---------------------------------------------------------------
    ADDITIONAL GLOBAL PARAMETERS
    --------------------------------------------------------------- */
/* Inflation / substitution analysis depth (Sections 15, 19). */
INFLATION_DEPTH : 6 $
/* Heat-trace / diffusion kernel time samples (Section 17). */
HEAT_TIMES : [0.01, 0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 50.0] $
/* Zeta function: number of terms in Euler-product approximation (Sec 16). */
ZETA_TERMS : 40 $
/* Acceptance-window sampling resolution (Section 18). */
WINDOW_SAMPLES : 20 $

/* Validate range parameters. */
if not integerp(N_GON_MIN)  or N_GON_MIN  < 3
    then error("N_GON_MIN must be an integer >= 3; got",  N_GON_MIN)$
if not integerp(N_GON_MAX)  or N_GON_MAX  < N_GON_MIN
    then error("N_GON_MAX must be >= N_GON_MIN; got",     N_GON_MAX)$
if not integerp(N_GON_STEP) or N_GON_STEP < 1
    then error("N_GON_STEP must be a positive integer; got", N_GON_STEP)$

print("------------------------------------------------------------")$
print("RANGE PARAMETERIZATION:")$
print("  N_GON_MIN  =", N_GON_MIN,
       "  N_GON_MAX =", N_GON_MAX,
       "  N_GON_STEP =", N_GON_STEP)$
print("------------------------------------------------------------")$

/* Accumulator for the per-n summary table (populated inside the loop). */
ngon_summary_table : [] $
/* Accumulator for cross-section correlation data (Section 20). */
cross_section_data : [] $

/* ---------------------------------------------------------------------
   SECTION 1.  Golden ratio phi and Q(sqrt(5))    [pentagon-specific]
   ---------------------------------------------------------------------

   These identities are pentagon-specific; we keep them as exact
   symbolic checks so any future code that imports the file can rely
   on phi, psi, and their identities.  Generalizations to other n
   involve different number fields:
        n=3,4,6 : Q              (degree 1)
        n=5,8,10,12 : real quadratic field over Q
        n=7,9,11,13,14 : higher-degree real cyclotomic fields
   --------------------------------------------------------------------- */

phi : (1 + sqrt(5))/2 $
psi : (1 - sqrt(5))/2 $   /* Galois conjugate of phi */

P1 : ratsimp(phi^2 - (phi + 1)) $
print("[Pentagon] P1: phi^2 - (phi+1) =", P1)$
if P1 # 0 then error("Identity P1 (phi^2 = phi+1) failed")$

P2 : ratsimp(1/phi - (phi - 1)) $
print("[Pentagon] P2: 1/phi - (phi-1) =", P2)$
if P2 # 0 then error("Identity P2 failed")$

P4 : ratsimp(phi*(phi - 1) - 1) $
print("[Pentagon] P4: phi*(phi-1) - 1 =", P4)$
if P4 # 0 then error("Identity P4 failed")$

P6a : ratsimp(phi + psi - 1) $
P6b : ratsimp(phi*psi - (-1)) $
print("[Pentagon] P6a: phi + psi - 1 =", P6a,
      "  (Vieta: sum of roots of x^2-x-1)")$
print("[Pentagon] P6b: phi*psi + 1 =", P6b,
      "  (Vieta: product of roots)")$
if P6a # 0 then error("Identity P6a failed")$
if P6b # 0 then error("Identity P6b failed")$

P7 : ratsimp(phi^2 + psi^2 - 3) $
print("[Pentagon] P7: phi^2 + psi^2 - 3 =", P7,
      "  (trace of phi^2 in Q(sqrt 5))")$
if P7 # 0 then error("Identity P7 failed")$

minpoly_phi(x) := x^2 - x - 1 $
print("[Pentagon] Minimal polynomial of phi:", minpoly_phi(x),
      "  value at phi:", ratsimp(minpoly_phi(phi)),
      ",  value at psi:", ratsimp(minpoly_phi(psi)))$
if ratsimp(minpoly_phi(phi)) # 0
   then error("Minimal polynomial check failed for phi")$
if ratsimp(minpoly_phi(psi)) # 0
   then error("Minimal polynomial check failed for psi")$

qsqrt5(a, b) := a + b*sqrt(5) $
qconj(z) := subst(-sqrt(5), sqrt(5), z) $
qnorm(z) := ratsimp(z * qconj(z)) $
qtrace(z) := ratsimp(z + qconj(z)) $

print("[Pentagon] qnorm(phi) =", qnorm(phi), " (expected -1)")$
if qnorm(phi) # -1 then error("Norm of phi expected to be -1")$
print("[Pentagon] qtrace(phi) =", qtrace(phi), " (expected 1)")$
if qtrace(phi) # 1 then error("Trace of phi expected to be 1")$

print("[Pentagon] Fibonacci-phi identity check (phi^n = F_n phi + F_{n-1}):")$
for nn:2 thru 8 do block(
    [lhs, rhs, diff_val],
    lhs : expand(phi^nn),
    rhs : fib(nn)*phi + fib(nn-1),
    diff_val : ratsimp(lhs - rhs),
    print("  n =", nn, ":  diff =", diff_val),
    if diff_val # 0 then error("Fibonacci-phi identity failed at n =", nn))$

print("[Pentagon] Binet formula check F_n = (phi^n - psi^n)/sqrt(5):")$
for nn:1 thru 8 do block(
    [binet, diff_val],
    binet : ratsimp((phi^nn - psi^nn)/sqrt(5)),
    diff_val : ratsimp(binet - fib(nn)),
    print("  n =", nn, ": F_n =", fib(nn), ", Binet =", binet,
          ", diff =", diff_val),
    if diff_val # 0 then error("Binet identity failed at n =", nn))$

/* ===================================================================
    HELPER FUNCTIONS (defined once, used inside the loop)
    =================================================================== */

/* General formulas as functions of n. */
interior_angle_ngon(n) := (n - 2)*%pi / n $
exterior_angle_ngon(n) := 2*%pi / n $
sum_interior_ngon(n)   := (n - 2)*%pi $

/* Rotation matrix by 2*pi/n. */
Rn_general(n) := matrix(
     [cos(2*%pi/n), -sin(2*%pi/n)],
     [sin(2*%pi/n),  cos(2*%pi/n)]) $

/* Vertex angle for vertex k of a regular n-gon (vertex 0 at top). */
vertex_angle_ngon(n, k) := %pi/2 + 2*%pi*k/n $

/* Holonomy helpers. */
holonomy_Zn(transitions, n) :=
     mod(lreduce("+", transitions, 0), n) $

order_Zn(g, n) := block([gg:mod(g, n), result],
      if gg = 0 then return(1),
      result : catch(
          for k:1 thru n do
              if mod(k*gg, n) = 0 then throw(k),
          n),
      result) $

/* Alexander-Orbach spectral dimension. */
d_spec_AO(d_eff, d_w) := 2 * d_eff / d_w $

/* Cyclic shift matrix. */
cyclic_shift_n(n) := genmatrix(
     lambda([i, j], if mod(j - i, n) = 1 then 1 else 0),
     n, n) $

/* CA rule-count formulas. */
num_outer_tot_rules(k)     := 2^(2*(k+1)) $
num_tot_rules(k)           := 2^(k+2) $
num_general_rules(k)       := 2^(2^(k+1)) $
num_outer_tot_rules_q(k,q) := q^(q * (k*(q-1) + 1)) $

/* Polygon geometry (unit side). */
circumradius_ngon(n) := 1/(2*sin(%pi/n)) $
apothem_ngon(n)      := 1/(2*tan(%pi/n)) $
area_ngon(n)         := (n/4) * cot(%pi/n) $
diagonal_ngon(n, k)  := sin(k*%pi/n)/sin(%pi/n) $

/* Acceptance-window helper. */
in_acceptance_window(pperp, rho) :=
     is(pperp . transpose(pperp) <= rho^2) $

/* Dimensional-flow interpolator. */
d_spec_flow(d_uv, d_ir, t, t0) :=
     d_ir + (d_uv - d_ir) * exp(-t/t0) $

/* n-independent Alexander-Orbach checks (run once before the loop). */
print("------------------------------------------------------------")$
print("Emergent-dimension framework (n-independent, run once):")$

cond : ratsimp(d_spec_AO(d_eff, 2) - d_eff) $
print("  d_spec(d_eff, d_w=2) - d_eff =", cond,
       "  (d_w=2 saturates Alexander-Orbach)")$
if cond # 0 then error("Alexander-Orbach saturation failed")$

dspec_deriv_dw   : diff(d_spec_AO(d_eff, d_w), d_w) $
dspec_deriv_deff : diff(d_spec_AO(d_eff, d_w), d_eff) $
print("  d/d(d_w)   d_spec =", dspec_deriv_dw)$
print("  d/d(d_eff) d_spec =", dspec_deriv_deff)$

msd_exponent(d_w)    := 2 / d_w $
p0_exponent(d_spec)  := d_spec / 2 $
print("  MSD exponent  2/d_w  at d_w=2.5 =", float(msd_exponent(2.5)))$
print("  P_0(t) decay exponent d_spec/2 at d_spec=1.5 =",
       float(p0_exponent(1.5)))$

test_dims(d_eff_val, d_w_val) :=
     block([d_s],
         d_s : float(d_spec_AO(d_eff_val, d_w_val)),
         print("  d_eff=", d_eff_val, " d_w=", d_w_val,
               " => d_spec=", d_s,
               if d_s < d_eff_val then "  (sub-diffusive)"
               else "  (not sub-diffusive)"),
         d_s) $

print("Dimension sanity table:")$
test_dims(2.5, 2.0) $
test_dims(2.5, 2.5) $
test_dims(2.5, 3.0) $
test_dims(2.8, 3.2) $
test_dims(2.2, 2.4) $

sier_deff  : log(3)/log(2) $
sier_dw    : log(5)/log(2) $
sier_dspec : float(d_spec_AO(sier_deff, sier_dw)) $
print("Sierpinski-gasket reference: d_eff=",
       float(sier_deff), " d_w=", float(sier_dw),
       " d_spec=", sier_dspec)$

in_paper_regime(d_eff_val, d_w_val) :=
     is(d_eff_val > 2) and is(d_eff_val < 3)
     and is(d_w_val > 2)
     and is(float(d_spec_AO(d_eff_val, d_w_val)) < d_eff_val) $

print("in_paper_regime(2.5, 2.4) =", in_paper_regime(2.5, 2.4),
       " (true expected)")$
print("in_paper_regime(2.0, 2.0) =", in_paper_regime(2.0, 2.0))$
print("in_paper_regime(2.5, 1.9) =", in_paper_regime(2.5, 1.9))$
print("in_paper_regime(3.1, 2.5) =", in_paper_regime(3.1, 2.5))$

print("Dimensional flow d_spec(t) interpolator (UV=2.5, IR=1.7):")$
for tt in [0.1, 1, 10, 100] do
     print("  t =", tt,
           ": d_spec(t) =", float(d_spec_flow(2.5, 1.7, tt, 10))) $
/* ---------------------------------------------------------------
    EXTENDED ALEXANDER-ORBACH ANALYSIS
    --------------------------------------------------------------- */
print("------------------------------------------------------------")$
print("Extended Alexander-Orbach / spectral-dimension analysis:")$
/* Full (d_eff, d_w) phase diagram on a grid. */
print("  Phase diagram: d_spec = 2*d_eff/d_w on grid")$
print("  (rows = d_eff 1.5..3.5 step 0.5; cols = d_w 1.5..4.0 step 0.5)")$
block([deff_vals, dw_vals, header, sep],
     deff_vals : [1.5, 2.0, 2.5, 3.0, 3.5],
     dw_vals   : [1.5, 2.0, 2.5, 3.0, 3.5, 4.0],
     header : pad_right("d_eff\\d_w", 12),
     for dw in dw_vals do header : sconcat(header, " | ", pad_right(string(dw), 7)),
     print("  ", header),
     sep : "  ------------",
     for dw in dw_vals do sep : sconcat(sep, "-+---------"),
     print(sep),
     for deff in deff_vals do block(
         [row],
         row : pad_right(string(deff), 12),
         for dw in dw_vals do block(
             [ds],
             ds  : float(d_spec_AO(deff, dw)),
             row : sconcat(row, " | ", pad_right(string(ds), 7))),
         print("  ", row))) $
/* Anomalous diffusion classification. */
classify_diffusion(d_w) :=
     if d_w < 2 then "super-diffusive"
     else if d_w = 2 then "normal"
     else "sub-diffusive (anomalous)" $
print("  Diffusion classification by d_w:")$
for dw in [1.5, 2.0, 2.5, 3.0, 4.0] do
     print("    d_w =", dw, "->", classify_diffusion(dw)) $
/* Spectral dimension from return probability exponent. */
/* P_0(t) ~ t^{-d_s/2}  =>  log P_0 / log t = -d_s/2 */
print("  Return-probability log-log slope (d_s/2) for d_s in {1.0..2.5}:")$
for ds in [1.0, 1.2, 1.5, 1.7, 2.0, 2.5] do
     print("    d_s =", ds, "  slope =", -ds/2,
           "  P_0(t=100) =", float(100^(-ds/2))) $
/* Crossover scale: when does UV regime end? */
/* d_spec(t) crosses d_eff at t = t_cross where exp(-t/t0) = (d_eff-d_ir)/(d_uv-d_ir) */
t_crossover(d_uv, d_ir, d_eff_val, t0) :=
     if d_uv = d_ir then inf
     else -t0 * log((d_eff_val - d_ir)/(d_uv - d_ir)) $
print("  Crossover times (UV=2.5, IR=1.7, t0=10):")$
for deff in [1.8, 2.0, 2.2, 2.4] do
     print("    d_eff =", deff, "  t_cross =",
           float(t_crossover(2.5, 1.7, deff, 10))) $
/* Integrated spectral weight. */
/* integral_0^Lambda rho(lambda) d lambda = d_spec/2 * integral lambda^{d_s/2-1} */
dos_integral(lam_max, B, d_spec) :=
     B * lam_max^(d_spec/2) / (d_spec/2) $
print("  DOS integrated weight (B=1, d_s=1.7, Lambda=1) =",
       float(dos_integral(1.0, 1.0, 1.7))) $
print("  DOS integrated weight (B=1, d_s=2.0, Lambda=1) =",
       float(dos_integral(1.0, 1.0, 2.0))) $

/* Volume-growth and MSD models (n-independent, run once). */
print("------------------------------------------------------------")$
print("Volume-growth and MSD models (n-independent, run once):")$

synthetic_volume(r, c, d) := c * r^d $

sample_radii   : makelist(rr, rr, 1, 10) $
sample_volumes : map(lambda([rr], float(synthetic_volume(rr, 2.5, 2.5))),
                      sample_radii) $
print("  Sample r       :", sample_radii)$
print("  Sample V(r)    :", sample_volumes)$
log_slope : float((log(sample_volumes[10]) - log(sample_volumes[1]))
                   / (log(10) - log(1))) $
print("  Recovered slope:", log_slope, " (expected 2.5)")$
if abs(log_slope - 2.5) > 1.0e-10
    then error("Slope recovery failed")$

msd_model(t, D_alpha, d_w) := D_alpha * t^(2/d_w) $
print("  MSD model (D=1, d_w=2.5):")$
for tt in [1, 10, 100, 1000] do
     print("    t =", tt, ":  MSD =", float(msd_model(tt, 1, 2.5))) $

p0_model(t, A, d_spec) := A * t^(-d_spec/2) $
print("  P_0(t) model (A=1, d_spec=1.7):")$
for tt in [1, 10, 100, 1000] do
     print("    t =", tt, ":  P_0 =", float(p0_model(tt, 1, 1.7))) $

dos_model(lam, B, d_spec) := B * lam^(d_spec/2 - 1) $
print("  DOS model (B=1, d_spec=1.7):")$
for lam in [0.001, 0.01, 0.1, 1] do
     print("    lambda =", lam, ":  rho =", float(dos_model(lam, 1, 1.7))) $
/* ---------------------------------------------------------------
    EXTENDED VOLUME-GROWTH ANALYSIS
    --------------------------------------------------------------- */
print("------------------------------------------------------------")$
print("Extended volume-growth analysis (n-independent):")$
/* Multi-scale volume growth: compare d_eff = 2.0, 2.5, 3.0 */
print("  V(r) = c * r^d_eff for d_eff in {2.0, 2.5, 3.0}, c=1:")$
for deff in [2.0, 2.5, 3.0] do block(
     [vols],
     vols : map(lambda([rr], float(synthetic_volume(rr, 1.0, deff))),
                [1, 2, 5, 10, 20, 50, 100]),
     print("    d_eff =", deff, ": V(r) =", vols)) $
/* Ratio V(2r)/V(r) = 2^d_eff (self-similarity). */
print("  Self-similarity ratio V(2r)/V(r) = 2^d_eff:")$
for deff in [1.5, 2.0, 2.5, 3.0] do
     print("    d_eff =", deff, "  ratio =", float(2^deff)) $
/* Effective dimension from two-point volume ratio. */
d_eff_from_ratio(V1, V2, r1, r2) :=
     log(V2/V1) / log(r2/r1) $
print("  Recovered d_eff from V(5)/V(1) with true d_eff=2.5:",
       float(d_eff_from_ratio(
           synthetic_volume(1,  1.0, 2.5),
           synthetic_volume(5,  1.0, 2.5), 1, 5))) $
print("  Recovered d_eff from V(10)/V(2) with true d_eff=2.5:",
       float(d_eff_from_ratio(
           synthetic_volume(2,  1.0, 2.5),
           synthetic_volume(10, 1.0, 2.5), 2, 10))) $
/* Lacunarity: variance of V(r) over random sub-windows. */
/* Lacunarity Lambda = (sigma_V / mu_V)^2; for self-similar sets Lambda ~ r^{-d_eff} */
lacunarity_model(r, d_eff_val) := r^(-d_eff_val) $
print("  Lacunarity model Lambda(r) = r^{-d_eff} (d_eff=2.5):")$
for rr in [1, 2, 5, 10, 20] do
     print("    r =", rr, "  Lambda =", float(lacunarity_model(rr, 2.5))) $
/* Multifractal spectrum: f(alpha) = d_eff for monofractal. */
/* For a monofractal, the Legendre transform of tau(q) = (q-1)*d_eff gives */
/* f(alpha) = d_eff at alpha = d_eff, zero elsewhere. */
tau_q_monofractal(q, d_eff_val) := (q - 1) * d_eff_val $
alpha_q_monofractal(q, d_eff_val) := d_eff_val $   /* d tau/d q */
f_alpha_monofractal(alpha, d_eff_val) :=
     if alpha = d_eff_val then d_eff_val else 0 $
print("  Monofractal tau(q) = (q-1)*d_eff (d_eff=2.5):")$
for q in [-2, -1, 0, 1, 2, 3] do
     print("    q =", q, "  tau(q) =", float(tau_q_monofractal(q, 2.5))) $

/* Z_2 spinor checks (n-independent, run once). */
print("------------------------------------------------------------")$
print("Spinor / holonomy checks (Z_2, n-independent, run once):")$
print("[Spinor case G = Z_2]:")$
print("  single-loop holonomy =", holonomy_Zn([1], 2),
       " (expected 1, i.e. tau = -1)")$
if holonomy_Zn([1], 2) # 1 then error("Z2 holonomy failed")$
print("  double-loop holonomy =", holonomy_Zn([1,1], 2),
       " (expected 0)")$
if holonomy_Zn([1,1], 2) # 0 then error("Z2 double-loop holonomy failed")$
print("  spinor loop order   =", order_Zn(1, 2),
       " (expected 2 -- double cover)")$
if order_Zn(1, 2) # 2 then error("Z2 spinor order != 2")$
/* ---------------------------------------------------------------
    EXTENDED HOLONOMY / FLAT-BUNDLE ANALYSIS (n-independent)
    --------------------------------------------------------------- */
print("------------------------------------------------------------")$
print("Extended flat-bundle / holonomy analysis (n-independent):")$
/* Holonomy group structure for Z_m covers. */
print("  Holonomy groups Z_m for m = 2..12:")$
for m : 2 thru 12 do block(
     [gens, orders],
     gens   : makelist(gg, gg, 0, m-1),
     orders : map(lambda([gg], order_Zn(gg, m)), gens),
     print("    Z_", m, ": generators and orders =",
           map(list, gens, orders))) $
/* Flat U(1) holonomy: phase accumulated around a vortex. */
/* For a Z_n vortex, the holonomy phase is exp(2 pi i k / n). */
flat_U1_holonomy(k, n) := exp(2*%pi*%i*k/n) $
print("  Flat U(1) holonomy phases exp(2 pi i k/n) for n=5:")$
for k : 0 thru 4 do
     print("    k =", k, ": phase =",
           rectform(flat_U1_holonomy(k, 5)),
           " ~ (", float(realpart(flat_U1_holonomy(k,5))),
           ",", float(imagpart(flat_U1_holonomy(k,5))), ")") $
/* Aharonov-Bohm phase: flux Phi through vortex, holonomy = exp(2 pi i Phi/Phi_0). */
AB_phase(flux_ratio) := exp(2*%pi*%i*flux_ratio) $
print("  Aharonov-Bohm phases for flux/flux_0 in {0, 1/5, 2/5, 3/5, 4/5, 1}:")$
for fr in [0, 1/5, 2/5, 3/5, 4/5, 1] do
     print("    Phi/Phi_0 =", fr, ": phase =",
           float(realpart(AB_phase(fr))), "+",
           float(imagpart(AB_phase(fr))), "i") $
/* Monodromy matrix for a 2D representation of Z_n. */
monodromy_2D(n) := matrix(
     [cos(2*%pi/n), -sin(2*%pi/n)],
     [sin(2*%pi/n),  cos(2*%pi/n)]) $
print("  Monodromy matrices for n = 3, 4, 5, 6:")$
for n_m in [3, 4, 5, 6] do block(
     [M, Mk, found_k],
     M : monodromy_2D(n_m),
     print("    n =", n_m, ": M =",
           matrix([float(M[1,1]), float(M[1,2])],
                  [float(M[2,1]), float(M[2,2])])),
     /* Find smallest k with M^k = I */
     Mk : ident(2),
     found_k : -1,
     for k : 1 thru 2*n_m do block(
         Mk : Mk . M,
         if found_k < 0 and
            abs(float(trigsimp(Mk[1,1])) - 1) < 1e-10 and
            abs(float(trigsimp(Mk[1,2])))     < 1e-10 and
            abs(float(trigsimp(Mk[2,1])))     < 1e-10 and
            abs(float(trigsimp(Mk[2,2])) - 1) < 1e-10
         then found_k : k),
     print("      M^k = I first at k =", found_k, " (expected", n_m, ")")) $
/* Berry phase for a closed loop in parameter space. */
/* For a 2-level system with gap Delta, Berry phase = pi * winding_number. */
berry_phase(winding) := %pi * winding $
print("  Berry phases for winding numbers 0..4:")$
for w : 0 thru 4 do
     print("    w =", w, ": gamma_B =", float(berry_phase(w)),
           " = ", w, "* pi") $
/* Chern number from holonomy: C = (1/2pi) * integral of Berry curvature. */
/* For Z_n holonomy, C = k/n for the k-th sector. */
chern_number(k, n) := k/n $
print("  Chern numbers for Z_5 sectors k = 0..4:")$
for k : 0 thru 4 do
     print("    k =", k, ": C =", chern_number(k, 5),
           " =", float(chern_number(k, 5))) $

/* ===================================================================
    MAIN SWEEP LOOP  over n = N_GON_MIN .. N_GON_MAX step N_GON_STEP
    =================================================================== */

print("============================================================")$
print("BEGIN SWEEP: n =", N_GON_MIN, "to", N_GON_MAX,
       "step", N_GON_STEP)$
print("============================================================")$
/* ---------------------------------------------------------------
    HELPER FUNCTIONS FOR ENHANCED SECTIONS (defined before loop)
    --------------------------------------------------------------- */
/* --- Section 15: Inflation / substitution helpers --- */
/* Apply a substitution matrix M to a vector v (tile counts). */
apply_substitution(M, v) := M . v $
/* Iterate substitution k times. */
iterate_substitution(M, v, k) := block(
     [result : v],
     for i : 1 thru k do result : apply_substitution(M, result),
     result) $
/* Perron-Frobenius eigenvalue (largest real eigenvalue, numeric). */
pf_eigenvalue_numeric(M) := block(
     [eigs, real_eigs],
     eigs      : map(lambda([e], float(e)), eigenvalues(M)[1]),
     real_eigs : sublist(eigs, lambda([e], imagpart(e) = 0)),
     if real_eigs = [] then error("No real eigenvalues found")
     else lmax(map(realpart, real_eigs))) $
/* Growth rate from iterated substitution (ratio of total counts). */
substitution_growth_rate(M, v, k) := block(
     [v1, v2, s1, s2],
     v1 : iterate_substitution(M, v, k),
     v2 : iterate_substitution(M, v, k+1),
     s1 : lreduce("+", list_matrix_entries(v1), 0),
     s2 : lreduce("+", list_matrix_entries(v2), 0),
     if s1 = 0 then error("Zero count at step k")
     else float(s2/s1)) $
/* --- Section 16: Zeta / L-series helpers --- */
/* Partial Euler product for zeta_Q(sqrt5)(s) up to norm N_max. */
/* Primes in Z[phi]: p splits if (5/p)=1, inert if (5/p)=-1, ramifies if p=5. */
legendre5(p) :=
     if mod(p, 5) = 0 then 0
     else if mod(p^2, 5) = 1 then 1   /* p = +-1 mod 5 => splits */
     else -1 $                          /* p = +-2 mod 5 => inert  */
/* Euler factor at rational prime p for zeta_{Q(sqrt5)}(s). */
euler_factor_qsqrt5(p, s) :=
     if legendre5(p) = 1  then 1/(1 - p^(-s))^2
     else if legendre5(p) = -1 then 1/(1 - p^(-2*s))
     else 1/(1 - p^(-s)) $   /* ramified: p = 5 */
/* Partial Dedekind zeta product over primes up to P_max. */
partial_dedekind_zeta(s, P_max) := block(
     [prod : 1.0, p],
     p : 2,
     while p <= P_max do block(
         prod : prod * float(euler_factor_qsqrt5(p, s)),
         p    : next_prime(p)),
     prod) $
/* Dirichlet L-function for the Kronecker symbol (5/.) */
/* L(s, chi_5) = sum_{n=1}^{N} chi_5(n) / n^s  (partial sum) */
kronecker5(n) :=
     if mod(n, 5) = 0 then 0
     else if member(mod(n, 5), [1, 4]) then 1
     else -1 $
partial_L_chi5(s, N_terms) :=
     float(sum(kronecker5(nn) / nn^s, nn, 1, N_terms)) $
/* --- Section 17: Heat-trace helpers --- */
/* Heat kernel on Z (1D lattice): K_t(x) = exp(-t) * I_x(t) */
/* Approximated by Gaussian for large t: K_t(x) ~ (4 pi t)^{-1/2} exp(-x^2/(4t)) */
heat_kernel_1D_gauss(x, t) := (1/sqrt(4*%pi*t)) * exp(-x^2/(4*t)) $
/* Heat trace on a finite cycle C_n: Tr(e^{-t L}) = sum_k exp(-t lambda_k) */
/* where lambda_k = 2 - 2 cos(2 pi k / n). */
heat_trace_cycle(n, t) :=
     float(sum(exp(-t * (2 - 2*cos(2*%pi*k/n))), k, 0, n-1)) $
/* Heat trace on star K_{1,n}: eigenvalues 0, 1 (x n-1), n+1. */
heat_trace_star(n, t) :=
     float(1 + (n-1)*exp(-t) + exp(-t*(n+1))) $
/* Spectral zeta function: zeta_L(s) = sum_{lambda_k > 0} lambda_k^{-s} */
/* For cycle C_n: */
spectral_zeta_cycle(n, s) :=
     float(sum(
         if k = 0 then 0
         else (2 - 2*cos(2*%pi*k/n))^(-s),
         k, 0, n-1)) $
/* Effective spectral dimension from heat-trace slope: */
/* Tr(e^{-t L}) ~ t^{-d_s/2}  =>  d_s = -2 * d(log Tr)/d(log t) */
heat_trace_slope(n, t1, t2) :=
     -2 * (log(heat_trace_cycle(n, t2)) - log(heat_trace_cycle(n, t1)))
          / (log(t2) - log(t1)) $
/* --- Section 18: Acceptance-window helpers --- */
/* Regular n-gon acceptance window: vertices at angles 2 pi k/n, radius rho. */
window_vertex(n, k, rho) :=
     [rho * cos(2*%pi*k/n), rho * sin(2*%pi*k/n)] $
/* Check if a 2D point p_perp lies inside the regular n-gon window of radius rho. */
/* Uses the half-plane intersection test. */
in_ngon_window(p_perp, n, rho) := block(
      [inside : true, px : p_perp[1], py : p_perp[2]],
     for k : 0 thru n-1 do block(
         [theta_k : 2*%pi*k/n,
          nx, ny, d],
         nx : float(cos(theta_k + %pi/2)),
         ny : float(sin(theta_k + %pi/2)),
         d  : float(rho * cos(%pi/n)),
          if nx*px + ny*py > d then inside : false),
     inside) $
/* Area of regular n-gon with inradius rho. */
ngon_window_area(n, rho) := n * rho^2 * tan(%pi/n) $
/* Density of accepted points: rho_accept = A_window / A_unit_cell */
/* For the pentagonal lattice, A_unit_cell = phi^2 * sin(2 pi/5). */
acceptance_density(n, rho, A_cell) :=
     float(ngon_window_area(n, rho)) / float(A_cell) $
/* --- Section 19: Fibonacci word / substitution helpers --- */
/* Fibonacci word substitution: a -> ab, b -> a */
fib_word_substitute(w) := block(
     [result : ""],
     for i : 1 thru slength(w) do block(
         [c : substring(w, i, i+1)],
         result : sconcat(result,
                          if c = "a" then "ab" else "a")),
     result) $
/* Generate Fibonacci word of depth k starting from "a". */
fib_word(k) := block(
     [w : "a"],
     for i : 1 thru k do w : fib_word_substitute(w),
     w) $
/* Count occurrences of character c in string w. */
count_char(w, c) := block(
     [cnt : 0],
     for i : 1 thru slength(w) do
         if substring(w, i, i+1) = c then cnt : cnt + 1,
     cnt) $
/* Thue-Morse sequence: t(n) = parity of number of 1s in binary(n). */
thue_morse(n) := mod(sum(
     block([bit : floor(n / 2^k)],
           mod(bit, 2)),
     k, 0, floor(log(max(n,1))/log(2))), 2) $
/* Rudin-Shapiro sequence: r(n) = (-1)^{number of 11 pairs in binary(n)}. */
rudin_shapiro(n) := block(
     [bits : [], m : n, cnt : 0],
     while m > 0 do (bits : cons(mod(m,2), bits), m : floor(m/2)),
     for i : 1 thru length(bits)-1 do
         if bits[i] = 1 and bits[i+1] = 1 then cnt : cnt + 1,
     (-1)^cnt) $

for N_GON : N_GON_MIN thru N_GON_MAX step N_GON_STEP do block(

   /* ---- local variables for this iteration ---- */
   [n_int, theta_int_n, theta_ext_n, theta_cen_n,
    cos_int_n, sin_int_n, cos_ext_n, sin_ext_n,
    k_flat_n, deficit_n, excess_n,
    k_loop_close_n, n_full_turns_n,
    Rn, detRn, trRn, Rn_power, Rn_power_floats, Rk,
    vertices_unit, vertices_unit_num,
    edges_unit, edges_unit_num, side_unit_R,
    fiber_order_n, h1, h2, h_comp, h_sum,
    sum_all_roots, sum_prim_roots_n, moebius_n, two_cos_n,
    Cn, In_id, charCn, Cn_power, Cn_orth,
    ones_n, eig_diag, Jn_mat, P_diag_n, P_perp_n, rank_perp_n,
    n_half, demo_perp, demo_norm2,
    R_unit_n, r_unit_n, A_unit_n,
    S0, conj_Rn, Rn_inv_p, conj_diff, Dn_order,
    A_star_n, D_star_n, L_star_n, ones_star,
    L_star_char, L_star_char_expected,
    A_Cn, D_Cn, L_Cn, Cn_eigs_num, Cn_gap,
    row_entry],

   /* ----------------------------------------------------------------
      Derived constants for this n
      ---------------------------------------------------------------- */
   n_int       : N_GON,
   theta_int_n : (n_int - 2)*%pi / n_int,
   theta_ext_n : 2*%pi / n_int,
   theta_cen_n : 2*%pi / n_int,
   cos_int_n   : cos(theta_int_n),
   sin_int_n   : sin(theta_int_n),
   cos_ext_n   : cos(theta_ext_n),
   sin_ext_n   : sin(theta_ext_n),

   k_flat_n  : floor(2*%pi / theta_int_n),
   deficit_n : ratsimp(2*%pi - k_flat_n * theta_int_n),
   excess_n  : ratsimp((k_flat_n + 1) * theta_int_n - 2*%pi),

   k_loop_close_n : 2*n_int / gcd(2*n_int, n_int - 2),
   n_full_turns_n : k_loop_close_n * (n_int - 2) / (2 * n_int),

   print("============================================================"),
   print("N_GON =", N_GON),
   print("============================================================"),

   print("Interior angle theta = (n-2)*pi/n =", theta_int_n,
         " =", float(theta_int_n*180/%pi), "deg"),
   print("Exterior angle = 2*pi/n          =", theta_ext_n,
         " =", float(theta_ext_n*180/%pi), "deg"),
   print("Flat-vertex count k_flat = floor(2*pi/theta) =", k_flat_n),
   print("Angular deficit at k_flat copies =", deficit_n,
         " =", float(deficit_n*180/%pi), "deg"),
   print("Angular excess at (k_flat+1) copies =", excess_n,
         " =", float(excess_n*180/%pi), "deg"),
   print("Smallest k>0 with k*theta = (integer)*2*pi:  k =",
         k_loop_close_n,
         ",  accumulating ", n_full_turns_n, " full turns"),

   /* ----------------------------------------------------------------
      SECTION 2.  Angles
      ---------------------------------------------------------------- */
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] Angle summary:"),
   print("  Interior angle  theta   =", theta_int_n,
         " =", float(theta_int_n), " rad =",
         float(theta_int_n*180/%pi), " deg"),
   print("  Exterior angle          =", theta_ext_n,
         " =", float(theta_ext_n), " rad"),
   print("  Central angle           =", theta_cen_n,
         " =", float(theta_cen_n), " rad"),
   print("  Sum of interior angles  =", sum_interior_ngon(n_int)),
   print("  Avg interior angle      =",
         ratsimp(sum_interior_ngon(n_int)/n_int)),

   if ratsimp(sum_interior_ngon(n_int)/n_int - theta_int_n) # 0
      then error("Average interior angle inconsistent at n =", n_int),

   print("[n = ", N_GON, "] k-copy vertex-loop accumulation:"),
   for kk:1 thru max(12, 2*N_GON) do block(
       [accum, n_turns, closes],
       accum   : ratsimp(kk * theta_int_n),
       n_turns : ratsimp(accum / (2*%pi)),
       closes  : if integerp(ratsimp(n_turns)) then "CLOSES" else "open",
       print("  k =", kk, ": sum =", accum, " = ",
             ratsimp(n_turns), "* 2pi  -- ", closes)),

   print("Smallest k with full closure: k =", k_loop_close_n,
         "  (number of full turns =", n_full_turns_n, ")"),

   /* ----------------------------------------------------------------
      SECTION 3.  Roots of unity
      ---------------------------------------------------------------- */
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] Cyclotomic / roots-of-unity structure:"),
   print("  x^n - 1 factorization =", factor(x^n_int - 1)),
   print("  totient(n) = phi(n)   =", totient(n_int),
         "  (degree of Q[zeta_n] over Q)"),
   print("  Real subfield degree  =",
         if n_int <= 2 then 1 else totient(n_int)/2,
         "  (degree of Q(zeta + zeta^{-1}) over Q)"),

   sum_all_roots : rectform(sum(exp(2*%pi*%i*k/n_int), k, 0, n_int - 1)),
   print("  Sum of all n-th roots (numeric) =", float(sum_all_roots),
         " (expected 0 for n >= 2)"),
   if abs(float(realpart(sum_all_roots))) > 1.0e-10
      or abs(float(imagpart(sum_all_roots))) > 1.0e-10
      then error("Sum of n-th roots not zero at n =", n_int),

   sum_prim_roots_n : rectform(
       sum(if gcd(k, n_int) = 1 then exp(2*%pi*%i*k/n_int) else 0,
           k, 1, n_int - 1)),
   moebius_n : moebius(n_int),
   print("  Sum of primitive n-th roots (numeric) =",
         float(sum_prim_roots_n),
         ",  Moebius mu(n) =", moebius_n),
   if abs(float(realpart(sum_prim_roots_n)) - moebius_n) > 1.0e-10
      or abs(float(imagpart(sum_prim_roots_n))) > 1.0e-10
      then error("Sum of primitive n-th roots != mu(n) at n =", n_int),

   two_cos_n : 2*cos(2*%pi/n_int),
   print("  2*cos(2*pi/n) =", two_cos_n, " =", float(two_cos_n),
         "  (algebraic integer in real cyclotomic subfield)"),

   if n_int = 5 then block(
       [check_72, check_144],
       check_72  : float(2*cos(2*%pi/5)  - ((sqrt(5)-1)/2)),
       check_144 : float(2*cos(4*%pi/5)  - (-(sqrt(5)+1)/2)),
       print("  [Pentagon] 2*cos(72)  - (sqrt 5 - 1)/2  ~", check_72),
       print("  [Pentagon] 2*cos(144) - (-(sqrt 5 + 1)/2) ~", check_144),
       if abs(check_72) > 1.0e-12 or abs(check_144) > 1.0e-12
          then error("Pentagon 2cos identities failed")),

   /* ----------------------------------------------------------------
      SECTION 4.  Rotation matrices and basis vectors
      ---------------------------------------------------------------- */
   Rn : Rn_general(n_int),
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] Rotation matrix R (by exterior angle 2*pi/n):"),
   print("  Symbolic:", Rn),
   print("  Numeric :",
         matrix([float(cos(2*%pi/n_int)), float(-sin(2*%pi/n_int))],
                [float(sin(2*%pi/n_int)), float( cos(2*%pi/n_int))])),

   detRn : ratsimp(trigsimp(determinant(Rn))),
   print("  det(R) =", detRn, " (expected 1)"),
   if ratsimp(detRn - 1) # 0 then error("det(R) != 1 at n =", n_int),

   trRn : ratsimp(mat_trace(Rn)),
   print("  tr(R) =", trRn, " = 2 cos(2*pi/n) =", float(trRn)),

   Rn_power : Rn,
   for kk:2 thru n_int do Rn_power : Rn_power . Rn,
   Rn_power_floats : matrix(
       [float(trigsimp(Rn_power[1,1])), float(trigsimp(Rn_power[1,2]))],
       [float(trigsimp(Rn_power[2,1])), float(trigsimp(Rn_power[2,2]))]),
   print("  R^n  numeric =", Rn_power_floats, " (expected I)"),
   for i:1 thru 2 do for j:1 thru 2 do
       if abs(Rn_power_floats[i,j] - if i=j then 1 else 0) > 1.0e-10
          then error("R^n != I at n =", n_int),

   print("[n = ", N_GON, "] R^k traces (should equal 2 cos(2 pi k/n)):"),
   Rk : ident(2),
   for kk:1 thru n_int do block(
       [trk, expected],
       Rk       : Rk . Rn,
       trk      : float(trigsimp(mat_trace(Rk))),
       expected : float(2*cos(2*%pi*kk/n_int)),
       print("  k =", kk, ": tr(R^k) =", trk, " ~ 2cos(",kk,
             "*2pi/n) =", expected),
       if abs(trk - expected) > 1.0e-10
          then error("R^k trace mismatch at k =", kk, "n =", n_int)),

   vertices_unit : makelist(
       [cos(vertex_angle_ngon(n_int, k)),
        sin(vertex_angle_ngon(n_int, k))],
       k, 0, n_int - 1),
   vertices_unit_num : makelist(
       [float(cos(vertex_angle_ngon(n_int, k))),
        float(sin(vertex_angle_ngon(n_int, k)))],
       k, 0, n_int - 1),

   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] PRIMARY BASIS VECTORS (unit circumradius,"),
   print("                            vertex 0 at top, CCW):"),
   for k:0 thru n_int - 1 do
       print("  v[", k, "] = ", vertices_unit[k+1],
             " ~ ", vertices_unit_num[k+1]),

   edges_unit : makelist(
       [vertices_unit[mod(k, n_int) + 1][1]
        - vertices_unit[mod(k - 1, n_int) + 1][1],
        vertices_unit[mod(k, n_int) + 1][2]
        - vertices_unit[mod(k - 1, n_int) + 1][2]],
       k, 1, n_int),
   edges_unit_num : makelist(
       [float(trigsimp(edges_unit[k][1])),
        float(trigsimp(edges_unit[k][2]))],
       k, 1, n_int),
   print("[n = ", N_GON, "] EDGE VECTORS e[k] = v[k] - v[k-1] (unit",
         "circumradius):"),
   for k:1 thru n_int do
       print("  e[", k, "] (numeric) = ", edges_unit_num[k]),

   side_unit_R : 2*sin(%pi/n_int),
   print("[n = ", N_GON, "] Side length (unit circumradius) s =",
         side_unit_R, " =", float(side_unit_R)),

   if n_int = 5 then block(
       [c72, s72, c36, s36, R72, R72p5],
       c72   : (sqrt(5) - 1)/4,
       s72   : sqrt(10 + 2*sqrt(5))/4,
       c36   : (sqrt(5) + 1)/4,
       s36   : sqrt(10 - 2*sqrt(5))/4,
       print("  [Pentagon] cos(72) - (phi-1)/2 =",
             ratsimp(c72 - (phi-1)/2)),
       print("  [Pentagon] cos(36) - phi/2     =",
             ratsimp(c36 - phi/2)),
       print("  [Pentagon] cos^2(72) + sin^2(72) =",
             ratsimp(c72^2 + s72^2)),
       print("  [Pentagon] cos^2(36) + sin^2(36) =",
             ratsimp(c36^2 + s36^2)),
       R72   : matrix([c72, -s72], [s72, c72]),
       R72p5 : ratsimp(R72 . R72 . R72 . R72 . R72),
       print("  [Pentagon] R72^5 (radcan'd) =",
             map(lambda([e], ratsimp(radcan(e))), R72p5)),
       print("  [Pentagon] R72^5 - I numeric =",
             map(lambda([e], float(e)), R72p5 - ident(2)))),

   /* ----------------------------------------------------------------
      SECTION 5.  Holonomy
      ---------------------------------------------------------------- */
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] Holonomy / spinor cover analysis:"),

   fiber_order_n : k_loop_close_n,
   print("  Natural single-vortex fiber order =", fiber_order_n),
   print("  Z_n single-vortex holonomy (n =", fiber_order_n, ") :",
         holonomy_Zn([1], fiber_order_n)),
   print("  Required loop count to close fiber :",
         order_Zn(1, fiber_order_n)),

   print("[Anyonic case G = Z_n with n =", n_int, "]:"),
   print("  single-vortex holonomy =", holonomy_Zn([1], n_int)),
   print("  loop order required    =", order_Zn(1, n_int),
         " (expected", n_int, ")"),
   if order_Zn(1, n_int) # n_int
      then error("Z_n anyonic order != n at n =", n_int),

   print("[n = ", N_GON, "] order(g, n) = n/gcd(g, n) check (g = 0..n-1):"),
   for gg:0 thru n_int - 1 do block(
       [o, expected],
       o        : order_Zn(gg, n_int),
       expected : if gg = 0 then 1 else n_int/gcd(gg, n_int),
       if o # expected
          then error("order(", gg, ",", n_int, ") =", o,
                     ", expected ", expected)),
   print("  All g in Z_n satisfy order(g) = n/gcd(g, n)."),

   h1     : holonomy_Zn([1, 2], n_int),
   h2     : holonomy_Zn([mod(3, n_int), mod(4, n_int)], n_int),
   h_comp : holonomy_Zn([1, 2, mod(3, n_int), mod(4, n_int)], n_int),
   h_sum  : mod(h1 + h2, n_int),
   print("[n = ", N_GON, "] Composite holonomy: Loop A =", h1,
         ", Loop B =", h2, ", Composite =", h_comp,
         ", Sum mod n =", h_sum),
   if h_comp # h_sum
      then error("Composite holonomy failed at n =", n_int),

   /* ----------------------------------------------------------------
      SECTION 7.  Cut-and-project lattice
      ---------------------------------------------------------------- */
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] Cut-and-project ambient lattice Z^n:"),

   Cn    : cyclic_shift_n(n_int),
   In_id : ident(n_int),
   print("  C_n cyclic shift matrix:"),
   print("  ", Cn),

   charCn : ratsimp(expand(charpoly(Cn, x))),
   print("  charpoly(C_n) =", charCn),
   if expand(charCn - (x^n_int - 1)) # 0
      and expand(charCn + (x^n_int - 1)) # 0
      and expand(charCn - (-1)^n_int*(x^n_int - 1)) # 0
      then error("char poly of C_n should be +/- (x^n - 1) at n =", n_int),

   Cn_power : Cn,
   for kk:2 thru n_int do Cn_power : Cn_power . Cn,
   if Cn_power # In_id then error("C_n^n != I_n at n =", n_int),
   print("  C_n^n = I_n verified."),

   Cn_orth : ratsimp(Cn . transpose(Cn)),
   if Cn_orth # In_id then error("C_n should be orthogonal at n =", n_int),
   print("  C_n is orthogonal: C_n . C_n^T = I_n."),

   ones_n   : genmatrix(lambda([i,j], 1), n_int, 1),
   eig_diag : ratsimp(Cn . ones_n - ones_n),
   print("  C_n . (1,...,1)^T - (1,...,1)^T =",
         transpose(eig_diag), " (should be zero)"),
   if eig_diag # genmatrix(lambda([i,j], 0), n_int, 1)
      then error("Diagonal eigenvector of C_n failed at n =", n_int),

   Jn_mat   : genmatrix(lambda([i, j], 1), n_int, n_int),
   P_diag_n : ratsimp((1/n_int) * Jn_mat),
   P_perp_n : ratsimp(In_id - P_diag_n),

   if ratsimp(P_diag_n . P_diag_n - P_diag_n) # zeromatrix(n_int, n_int)
      then error("P_diag not idempotent at n =", n_int),
   if ratsimp(P_perp_n . P_perp_n - P_perp_n) # zeromatrix(n_int, n_int)
      then error("P_perp not idempotent at n =", n_int),
   if ratsimp(P_diag_n . P_perp_n) # zeromatrix(n_int, n_int)
      then error("Projectors not orthogonal at n =", n_int),
   if ratsimp(P_diag_n + P_perp_n - In_id) # zeromatrix(n_int, n_int)
      then error("Projectors don't sum to I_n at n =", n_int),
   print("  Projectors P_diag, P_perp: idempotent, orthogonal, complete."),

   rank_perp_n : mat_trace(P_perp_n),
   print("  dim(orthogonal complement of diagonal) =", rank_perp_n,
         " = n - 1"),
   if rank_perp_n # n_int - 1
      then error("Perp projector rank should be n - 1 at n =", n_int),

   print("[n = ", N_GON, "] REAL-INVARIANT 2D SUBSPACES V_k of C_n:"),
   n_half : floor((n_int - 1)/2),
   for k:1 thru n_half do block(
       [u_vec, w_vec, u_num, w_num, theta_k],
       theta_k : 2*%pi*k/n_int,
       u_vec   : makelist(cos(2*%pi*k*j/n_int), j, 0, n_int - 1),
       w_vec   : makelist(sin(2*%pi*k*j/n_int), j, 0, n_int - 1),
        u_num   : map(lambda([xx], float(xx)), u_vec),
        w_num   : map(lambda([xx], float(xx)), w_vec),
       print("  V_", k, " (rotation angle 2 pi *", k, "/n =",
             float(theta_k), "):"),
       print("    u_", k, " (cos) =", u_num),
       print("    w_", k, " (sin) =", w_num)),

   if mod(n_int, 2) = 0 then block(
       [alt_vec],
       alt_vec : makelist((-1)^j, j, 0, n_int - 1),
       print("  alternating-sign eigenvector (eigenvalue -1) =", alt_vec)),

   print("  Physical 2D plane E_|| recommended = V_1."),
   print("  Perp subspace E_perp = (diagonal) +",
         n_half - 1, "* 2D pieces +",
         if mod(n_int, 2) = 0 then "1 alternating direction"
         else "0 alternating directions",
         "  =>  dim(E_perp) =", n_int - 2),

   demo_perp : transpose(genmatrix(
       lambda([i, j],
              if i = 1 then 1 - 1/n_int else -1/n_int),
       n_int, 1)),
    demo_norm2 : ratsimp(sum(demo_perp[1, j]^2, j, 1, n_int)),
   print("  ||pi_perp(e_1)||^2 =", demo_norm2, " (= (n-1)/n)"),
   if ratsimp(demo_norm2 - (n_int - 1)/n_int) # 0
      then error("perp-projection norm wrong at n =", n_int),

   /* ----------------------------------------------------------------
      SECTION 8.  CA rule counts
      ---------------------------------------------------------------- */
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] CA rule-count analysis (degree k = n):"),
   print("  Binary outer-totalistic n-neighbor rules =",
         num_outer_tot_rules(n_int),
         " (= 2^", 2*(n_int+1), ")"),
   print("  Binary fully totalistic n-neighbor rules =",
         num_tot_rules(n_int),
         " (= 2^", n_int+2, ")"),
   print("  Binary fully general    n-neighbor rules =",
         num_general_rules(n_int),
         " (= 2^(2^", n_int+1, "))"),

   if num_outer_tot_rules_q(n_int, 2) # num_outer_tot_rules(n_int)
      then error("Generic OT formula disagrees with binary case at n =",
                 n_int),

   if n_int = 5 then block(
       if num_outer_tot_rules(5) # 4096
          then error("Rule count for 5-regular OT family failed"),
       print("  [Pentagon] 4096 outer-totalistic rules confirmed.")),

   print("  Reference: 4-OT =", num_outer_tot_rules(4),
         ", 8-OT =", num_outer_tot_rules(8)),

   /* ----------------------------------------------------------------
      SECTION 9.  Polygon geometry
      ---------------------------------------------------------------- */
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] Polygon geometry (side s = 1):"),

   R_unit_n : circumradius_ngon(n_int),
   r_unit_n : apothem_ngon(n_int),
   A_unit_n : area_ngon(n_int),

   print("  Circumradius R =", R_unit_n, " ~", float(R_unit_n)),
   print("  Apothem      r =", r_unit_n, " ~", float(r_unit_n)),
   print("  R/r            =", ratsimp(R_unit_n/r_unit_n),
         " = sec(pi/n) =", float(sec(%pi/n_int))),
   print("  Area A         =", A_unit_n, " ~", float(A_unit_n)),

   print("  Diagonals (k-step chord lengths, unit side):"),
   for k:1 thru n_int - 1 do
       print("    d[", k, "] =", ratsimp(diagonal_ngon(n_int, k)),
             " ~", float(diagonal_ngon(n_int, k))),

   if n_int = 5 then block(
       [d2_sq, phi_sq, diag_check],
       d2_sq      : ratsimp((3 + sqrt(5))/2),
       phi_sq     : ratsimp(phi^2),
       diag_check : ratsimp(d2_sq - phi_sq),
       if diag_check # 0
          then error("Pentagon diagonal/side = phi failed"),
       print("  [Pentagon] (d/s)^2 = phi^2 confirmed")),

   if n_int = 5 then
       if ratsimp(1/phi^2 - (2 - phi)) # 0
          then error("1/phi^2 = 2 - phi failed")
       else print("  [Pentagon] pentagram self-similarity: 1/phi^2 = 2 - phi"),

   /* ----------------------------------------------------------------
      SECTION 11.  Dihedral group D_n
      ---------------------------------------------------------------- */
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] Dihedral group D_n symmetries:"),

   S0 : matrix([1, 0], [0, -1]),
   if ratsimp(S0 . S0) # ident(2)
      then error("Reflection should be involution at n =", n_int),
   if determinant(S0) # -1
      then error("Reflection determinant should be -1 at n =", n_int),

   conj_Rn  : ratsimp(trigsimp(S0 . Rn . S0)),
   Rn_inv_p : ratsimp(trigsimp(invert(Rn))),
   conj_diff : ratsimp(trigsimp(conj_Rn - Rn_inv_p)),
   print("  S0 . R . S0 - R^{-1} (symbolic):", conj_diff),
   for i:1 thru 2 do for j:1 thru 2 do
       if abs(float(trigsimp(conj_diff[i,j]))) > 1.0e-10
          then error("Dihedral relation S R S = R^{-1} failed at n =", n_int),
   print("  S R S = R^{-1} verified."),

   Dn_order : 2 * n_int,
   print("  |D_n| =", Dn_order),

   if mod(n_int, 2) = 1 then block(
       [num_classes],
       num_classes : 1 + (n_int - 1)/2 + 1,
       print("  Conjugacy classes (n odd): 1 + (n-1)/2 + 1 =",
             num_classes,
             " (sum of sizes 1 + (n-1)/2 * 2 + n =",
             1 + (n_int - 1) + n_int, "= 2n)")),
   if mod(n_int, 2) = 0 then block(
       [num_classes],
       num_classes : 1 + (n_int - 2)/2 + 1 + 2,
       print("  Conjugacy classes (n even): 1 + (n-2)/2 + 1 + 2 =",
             num_classes,
             " (sum of sizes 1 + (n-2)/2 * 2 + 1 + 2*(n/2) =",
             1 + (n_int - 2) + 1 + n_int, "= 2n)")),

   if mod(n_int, 2) = 1 then block(
       [sumsq],
       sumsq : 2*1 + (n_int - 1)/2 * 4,
       print("  Irrep dim^2 sum (n odd): 2 + 2(n-1) =", sumsq,
             "(should = 2n =", 2*n_int, ")"),
       if sumsq # 2*n_int then error("D_n irrep dim sum failed at n =", n_int)),
   if mod(n_int, 2) = 0 then block(
       [sumsq],
       sumsq : 4*1 + (n_int - 2)/2 * 4,
       print("  Irrep dim^2 sum (n even): 4 + 2(n-2) =", sumsq,
             "(should = 2n =", 2*n_int, ")"),
       if sumsq # 2*n_int then error("D_n irrep dim sum failed at n =", n_int)),

   /* ----------------------------------------------------------------
      SECTION 13.  Graph Laplacians
      ---------------------------------------------------------------- */
   print("------------------------------------------------------------"),
   print("[n = ", N_GON, "] Graph Laplacians (star and cycle):"),

   A_star_n : genmatrix(
       lambda([i, j],
              if (i = 1 and j > 1) or (j = 1 and i > 1) then 1 else 0),
       n_int + 1, n_int + 1),
   D_star_n : genmatrix(
       lambda([i, j],
              if i = j then (if i = 1 then n_int else 1) else 0),
       n_int + 1, n_int + 1),
   L_star_n : D_star_n - A_star_n,

   ones_star : genmatrix(lambda([i,j], 1), n_int + 1, 1),
   if ratsimp(L_star_n . ones_star) # zeromatrix(n_int + 1, 1)
      then error("Star Laplacian should annihilate (1,...,1)^T at n =", n_int),

    L_star_char_expected : lam * (lam - 1)^(n_int - 1) * (lam - (n_int + 1)),
   L_star_char : expand(charpoly(L_star_n, lam)),
    if expand(L_star_char - expand(L_star_char_expected)) # 0
       and expand(L_star_char + expand(L_star_char_expected)) # 0
       then block(
           print("  ACTUAL  charpoly =", L_star_char),
           print("  EXPECT  charpoly =", expand(L_star_char_expected)),
           error("Star Laplacian spectrum mismatch at n =", n_int)),
   print("  Star K_{1,", n_int, "}: spectrum = {0,",
          "1 (multiplicity", n_int - 1, "),", n_int + 1, "}"),
   print("  Spectral gap = 1 (positive => connected)."),

   A_Cn : genmatrix(
       lambda([i, j],
              if mod(j - i, n_int) = 1 or mod(i - j, n_int) = 1
                then 1 else 0),
       n_int, n_int),
   D_Cn : 2 * ident(n_int),
   L_Cn : D_Cn - A_Cn,

   print("  Cycle C_", n_int, " Laplacian eigenvalues = 2 - 2 cos(2 pi k/n):"),
   Cn_eigs_num : makelist(float(2 - 2*cos(2*%pi*k/n_int)), k, 0, n_int - 1),
   print("  ", Cn_eigs_num),

   Cn_gap : float(2 - 2*cos(2*%pi/n_int)),
   print("  Spectral gap of C_", n_int, " =", Cn_gap),

   if n_int = 5 then block(
       [Lcp],
       Lcp : ratsimp(expand(charpoly(L_Cn, lam))),
       print("  [Pentagon] charpoly(L_{C_5}) =", Lcp)),

   /* ----------------------------------------------------------------
      Accumulate summary row for this n
      ---------------------------------------------------------------- */
   row_entry : [
       n_int,
       float(theta_int_n * 180 / %pi),   /* interior angle (deg)      */
       k_flat_n,                          /* flat-vertex count         */
       float(deficit_n * 180 / %pi),      /* deficit (deg)             */
       k_loop_close_n,                    /* loop-closure index        */
       n_full_turns_n,                    /* full turns at closure     */
       totient(n_int),                    /* phi(n)                    */
       float(area_ngon(n_int)),           /* area (unit side)          */
       float(circumradius_ngon(n_int)),   /* circumradius (unit side)  */
       num_outer_tot_rules(n_int),        /* OT rule count             */
       float(Cn_gap)                      /* C_n spectral gap          */
   ],
   ngon_summary_table : endcons(row_entry, ngon_summary_table),

    print("  [n =", N_GON, "] iteration complete."),

    /* ----------------------------------------------------------------
       SECTION 15.  Inflation / substitution system
       ---------------------------------------------------------------- */
    print("------------------------------------------------------------"),
    print("[n = ", N_GON, "] Inflation / substitution system:"),


    block(
        [M_sub, v0, inflation_ratio_n, pf_eig, growth_rates],
         /* For a regular n-gon tiling, the substitution matrix encodes how
            each tile type inflates into smaller tiles.  For the pentagon (n=5)
            this is the Penrose substitution; for general n we use the
            "star-polygon" substitution where each n-gon inflates into
            phi_n copies of itself and (n-2) "kite" pieces.
            Here we use a simplified 2-type model: type A (n-gon) and
            type B (complementary rhombus / kite). */
         /* Assign M_sub based on n_int BEFORE any use. */
         M_sub : if n_int = 3
                 then ident(2)
                 else if n_int = 5
                 then matrix([2, 1], [1, 1])
                 else matrix([2, 1], [n_int - 3, n_int - 4]),

        /* Substitution matrix for n-gon + rhombus system.
          For n=5 we use the Penrose (kite/dart) substitution matrix
          M = [[2,1],[1,1]] whose Perron-Frobenius eigenvalue is phi^2
          (charpoly x^2 - 3x + 1, roots phi^2 and 1/phi^2).
          For general n >= 4 we use M = [[2, 1], [n-3, n-4]] which
          reduces to the Penrose case at n=5.
          For n=3: identity (equilateral triangle tiles flat). */


        v0 : matrix([1], [0]),   /* start with one A tile */

        print("  Substitution matrix M =", M_sub),
        print("  Initial vector v0 =", transpose(v0)),

        /* Iterate substitution */
        for depth : 1 thru INFLATION_DEPTH do block(
            [vk],
            vk : iterate_substitution(M_sub, v0, depth),
            print("  depth =", depth, ": tile counts =",
                  transpose(vk),
                  "  total =",
                  lreduce("+", list_matrix_entries(vk), 0))),

        /* Perron-Frobenius eigenvalue = inflation ratio */
        if n_int = 3
        then inflation_ratio_n : 1
        else inflation_ratio_n : pf_eigenvalue_numeric(M_sub),
        print("  Perron-Frobenius eigenvalue (inflation ratio) =",
              inflation_ratio_n),

        /* For n=5, inflation ratio should be phi^2 */
        if n_int = 5 then block(
            [phi2_num : float(phi^2)],
            print("  [Pentagon] phi^2 =", phi2_num,
                  "  PF eigenvalue =", inflation_ratio_n,
                  "  diff =", abs(inflation_ratio_n - phi2_num)),
            if abs(inflation_ratio_n - phi2_num) > 1.0e-8
               then error("Pentagon inflation ratio != phi^2")),

        /* Growth rates at successive depths */
        print("  Tile-count growth rates (ratio at successive depths):"),
        growth_rates : makelist(
            substitution_growth_rate(M_sub, v0, d),
            d, 1, INFLATION_DEPTH - 1),
        print("  ", growth_rates),
        print("  Converging to PF eigenvalue:", inflation_ratio_n),

        /* Fractal dimension from inflation: d_f = log(N_tiles) / log(inflation_ratio) */
        /* where N_tiles = total tiles after one inflation step starting from 1 tile. */
        if n_int # 3 then block(
            [N_tiles_1, d_frac],
            N_tiles_1 : lreduce("+",
                list_matrix_entries(apply_substitution(M_sub, v0)), 0),
            d_frac : if inflation_ratio_n > 1
                     then float(log(N_tiles_1) / log(inflation_ratio_n))
                     else 0,
            print("  Tiles after 1 inflation step:", N_tiles_1),
            print("  Fractal dimension estimate d_f = log(N)/log(lambda) =",
                  d_frac))),

    /* ----------------------------------------------------------------
       SECTION 16.  Zeta functions and L-series
       ---------------------------------------------------------------- */
    print("------------------------------------------------------------"),
    print("[n = ", N_GON, "] Zeta functions and L-series:"),

    block(
        [zeta_2_partial, zeta_3_partial, L_chi_2, L_chi_3,
         riemann_zeta_2, riemann_zeta_3],

        /* Partial Riemann zeta sums for reference */
        riemann_zeta_2 : float(sum(1/nn^2, nn, 1, ZETA_TERMS)),
        riemann_zeta_3 : float(sum(1/nn^3, nn, 1, ZETA_TERMS)),
        print("  Partial Riemann zeta(2) (", ZETA_TERMS, " terms) =",
              riemann_zeta_2, "  (exact pi^2/6 =", float(%pi^2/6), ")"),
        print("  Partial Riemann zeta(3) (", ZETA_TERMS, " terms) =",
              riemann_zeta_3, "  (Apery's constant ~ 1.202)"),

        /* Dedekind zeta for Q(sqrt(5)) at s=2 and s=3 */
        zeta_2_partial : partial_dedekind_zeta(2, 30),
        zeta_3_partial : partial_dedekind_zeta(3, 30),
        print("  Partial Dedekind zeta_{Q(sqrt5)}(2) (primes <= 30) =",
              zeta_2_partial),
        print("  Partial Dedekind zeta_{Q(sqrt5)}(3) (primes <= 30) =",
              zeta_3_partial),
        print("  Ratio zeta_{Q(sqrt5)}(2) / zeta(2)^2 ~ (class number factor):"),
        print("    =", zeta_2_partial / riemann_zeta_2^2),

        /* Dirichlet L-function L(s, chi_5) */
        L_chi_2 : partial_L_chi5(2, ZETA_TERMS),
        L_chi_3 : partial_L_chi5(3, ZETA_TERMS),
        print("  Partial L(2, chi_5) (", ZETA_TERMS, " terms) =", L_chi_2),
        print("  Partial L(3, chi_5) (", ZETA_TERMS, " terms) =", L_chi_3),
        print("  Exact L(2, chi_5) = pi^2 / (5 sqrt(5)) ~",
              float(%pi^2 / (5*sqrt(5)))),

        /* Functional equation check: xi(s) = xi(1-s) */
        /* For L(s, chi_5): completed function satisfies xi(s) = xi(1-s). */
        /* Numerically: L(s) / L(1-s) = (pi/5)^{s-1/2} * Gamma(1-s)/Gamma(s) * chi_factor */
        print("  Functional equation (numerical check at s=0.7):"),
        block(
            [s_val : 0.7,
             L_s, L_1ms, ratio],
            L_s   : partial_L_chi5(s_val,   ZETA_TERMS),
            L_1ms : partial_L_chi5(1-s_val, ZETA_TERMS),
            ratio : if abs(L_1ms) > 1e-15 then L_s / L_1ms else 0,
            print("    L(0.7) =", L_s, "  L(0.3) =", L_1ms,
                  "  ratio =", ratio)),

        /* Spectral zeta for cycle C_n */
        if n_int >= 3 then block(
            [sz2, sz3],
            sz2 : spectral_zeta_cycle(n_int, 2),
            sz3 : spectral_zeta_cycle(n_int, 3),
            print("  Spectral zeta of C_", n_int, " at s=2:", sz2),
            print("  Spectral zeta of C_", n_int, " at s=3:", sz3))),

    /* ----------------------------------------------------------------
       SECTION 17.  Diffusion kernel and heat trace
       ---------------------------------------------------------------- */
    print("------------------------------------------------------------"),
    print("[n = ", N_GON, "] Diffusion kernel and heat trace:"),

    block(
        [ht_cycle, ht_star, slopes],

        print("  Heat trace Tr(e^{-t L}) for cycle C_", n_int, ":"),
        ht_cycle : map(lambda([t], heat_trace_cycle(n_int, t)), HEAT_TIMES),
        for i : 1 thru length(HEAT_TIMES) do
            print("    t =", HEAT_TIMES[i], ":  Tr =", ht_cycle[i]),

        print("  Heat trace Tr(e^{-t L}) for star K_{1,", n_int, "}:"),
        ht_star : map(lambda([t], heat_trace_star(n_int, t)), HEAT_TIMES),
        for i : 1 thru length(HEAT_TIMES) do
            print("    t =", HEAT_TIMES[i], ":  Tr =", ht_star[i]),

        /* Effective spectral dimension from log-log slope of heat trace */
        print("  Effective d_s from heat-trace slope (cycle C_", n_int, "):"),
        slopes : makelist(
            heat_trace_slope(n_int, HEAT_TIMES[i], HEAT_TIMES[i+1]),
            i, 1, length(HEAT_TIMES)-1),
        for i : 1 thru length(slopes) do
            print("    t in [", HEAT_TIMES[i], ",", HEAT_TIMES[i+1], "]:  d_s ~",
                  slopes[i]),

        /* Short-time asymptotics: Tr ~ (4 pi t)^{-d/2} * Vol */
        /* For a 2D manifold: Tr ~ (4 pi t)^{-1} * Area + (1/6) * (4 pi t)^{-1/2} * Length + ... */
        print("  Short-time heat-trace asymptotics (2D Weyl law):"),
        print("    Tr(e^{-tL}) ~ (4 pi t)^{-1} * A + (4 pi t)^{-1/2} * B/6 + ..."),
        for t_val in [0.01, 0.05, 0.1] do block(
            [weyl_2d],
            weyl_2d : float(1/(4*%pi*t_val)),
            print("    t =", t_val, ":  Weyl leading term =", weyl_2d)),

        /* Heat kernel on the physical 2D plane (Gaussian approximation) */
        print("  Heat kernel K_t(x) on R^2 (Gaussian, d=2):"),
        for t_val in [0.1, 1.0, 10.0] do block(
            [K0, K1, K2],
            K0 : float((1/(4*%pi*t_val))),
            K1 : float((1/(4*%pi*t_val)) * exp(-1/(4*t_val))),
            K2 : float((1/(4*%pi*t_val)) * exp(-4/(4*t_val))),
            print("    t =", t_val, ":  K(0) =", K0,
                  "  K(1) =", K1, "  K(2) =", K2)),

        /* Spectral zeta from heat trace via Mellin transform */
        /* zeta_L(s) = (1/Gamma(s)) * integral_0^inf t^{s-1} (Tr e^{-tL} - dim_ker) dt */
        print("  Spectral zeta of C_", n_int, " at s=1,2,3:"),
        for s_val in [1, 2, 3] do
            print("    s =", s_val, ":  zeta_L =",
                  spectral_zeta_cycle(n_int, s_val))),

    /* ----------------------------------------------------------------
       SECTION 18.  Acceptance-window geometry
       ---------------------------------------------------------------- */
    print("------------------------------------------------------------"),
    print("[n = ", N_GON, "] Acceptance-window (cut-and-project) geometry:"),

    block(
        [rho_vals, A_cell_n, densities, sample_pts, accepted],

        /* Unit-cell area for the n-gon lattice (approximate). */
        /* For the pentagonal lattice: A_cell = phi^2 * sin(2 pi/5). */
        A_cell_n : float(area_ngon(n_int)),

        print("  Acceptance window: regular ", n_int, "-gon in perp space"),
        print("  Unit-cell area A_cell ~", A_cell_n),

        /* Acceptance density for various window radii */
        rho_vals : [0.5, 1.0, 1.5, 2.0],
        print("  Acceptance densities (rho = window inradius):"),
        densities : map(
            lambda([rho], acceptance_density(n_int, rho, A_cell_n)),
            rho_vals),
        for i : 1 thru length(rho_vals) do
            print("    rho =", rho_vals[i], ":  density =", densities[i]),

        /* Window vertices */
        print("  Window vertices (rho = 1):"),
        for k : 0 thru n_int - 1 do
            print("    v[", k, "] = (",
                  float(cos(2*%pi*k/n_int)), ",",
                  float(sin(2*%pi*k/n_int)), ")"),

        /* Sample points in perp space and test acceptance */
        print("  Sample acceptance test (rho = 1, grid [-1.5, 1.5]^2,",
              WINDOW_SAMPLES, "x", WINDOW_SAMPLES, " grid):"),
        block(
             [grid_step : 3.0 / WINDOW_SAMPLES, cnt_in : 0, cnt_tot : 0],
            for ix : 0 thru WINDOW_SAMPLES - 1 do
                for iy : 0 thru WINDOW_SAMPLES - 1 do block(
                     [gpx : -1.5 + ix * grid_step,
                      gpy : -1.5 + iy * grid_step],
                    cnt_tot : cnt_tot + 1,
                     if in_ngon_window([gpx, gpy], n_int, 1.0)
                       then cnt_in : cnt_in + 1),
            print("    Points inside window:", cnt_in, "/", cnt_tot,
                  "  fraction =", float(cnt_in / cnt_tot)),
            print("    Expected fraction (A_window/A_box) =",
                  float(ngon_window_area(n_int, 1.0) / 9.0))),

        /* Inflation of acceptance window: under inflation by phi, */
        /* the window scales by 1/phi (internal space contracts). */
        if n_int = 5 then block(
            [rho_inflated],
            rho_inflated : float(1.0 / phi),
            print("  [Pentagon] After one inflation step, window radius:",
                  rho_inflated, " (= 1/phi)"),
            print("  [Pentagon] Window area ratio (inflated/original) =",
                  float(ngon_window_area(5, rho_inflated) /
                        ngon_window_area(5, 1.0)),
                  " (expected 1/phi^2 =", float(1/phi^2), ")"))),

    /* ----------------------------------------------------------------
       SECTION 19.  Fibonacci word and substitution sequences
       ---------------------------------------------------------------- */
    print("------------------------------------------------------------"),
    print("[n = ", N_GON, "] Substitution sequences and symbolic dynamics:"),

    block(
        [fw, len_fw, cnt_a, cnt_b, ratio_ab,
         tm_seq, rs_seq],

        /* Fibonacci word */
        print("  Fibonacci word (substitution a->ab, b->a):"),
        for depth : 0 thru min(INFLATION_DEPTH, 8) do block(
            [w],
            w : fib_word(depth),
            print("    depth =", depth, ": length =", slength(w),
                  "  word =",
                  if slength(w) <= 30 then w
                  else sconcat(substring(w, 1, 28), "..."))),

        /* a/b ratio converges to phi */
        fw : fib_word(INFLATION_DEPTH),
        len_fw : slength(fw),
        cnt_a  : count_char(fw, "a"),
        cnt_b  : count_char(fw, "b"),
        ratio_ab : if cnt_b > 0 then float(cnt_a / cnt_b) else inf,
        print("  At depth", INFLATION_DEPTH, ": |a| =", cnt_a,
              "  |b| =", cnt_b,
              "  |a|/|b| =", ratio_ab,
              "  phi =", float(phi)),
        if abs(ratio_ab - float(phi)) > 0.1
           then print("  WARNING: a/b ratio not yet close to phi at depth",
                      INFLATION_DEPTH)
           else print("  a/b ratio converging to phi: OK"),

        /* Fibonacci word complexity: p(n) = n + 1 (Sturmian) */
        print("  Fibonacci word is Sturmian: complexity p(n) = n + 1"),
        print("  (minimum complexity for an aperiodic sequence)"),

        /* Thue-Morse sequence (relevant for n=4 square lattice) */
        print("  Thue-Morse sequence (first 16 terms):"),
        tm_seq : makelist(thue_morse(k), k, 0, 15),
        print("  ", tm_seq),

        /* Rudin-Shapiro sequence */
        print("  Rudin-Shapiro sequence (first 16 terms):"),
        rs_seq : makelist(rudin_shapiro(k), k, 0, 15),
        print("  ", rs_seq),

        /* For n=5: the Fibonacci word encodes the long/short edge sequence
           in the Penrose tiling along any line. */
        if n_int = 5 then block(
            [w6 : fib_word(6)],
            print("  [Pentagon] Fibonacci word depth 6 (Penrose edge sequence):"),
            print("    ", w6),
            print("    Length =", slength(w6),
                  "  (= F_8 =", fib(8), ")"),
            if slength(w6) # fib(8)
               then error("Fibonacci word length != F_8 at depth 6")),

        /* Substitution matrix for Fibonacci word: [[1,1],[1,0]] */
        block(
            [M_fib, eigs_fib],
            M_fib : matrix([1, 1], [1, 0]),
            print("  Fibonacci substitution matrix:", M_fib),
            print("  Characteristic polynomial:",
                  expand(charpoly(M_fib, x))),
            if expand(charpoly(M_fib, x) - (x^2 - x - 1)) # 0
               then error("Fibonacci substitution charpoly wrong"),
            print("  Eigenvalues: phi =", float(phi),
                  " and psi =", float(psi),
                  " (roots of x^2 - x - 1)"),
            print("  PF eigenvalue =", pf_eigenvalue_numeric(M_fib),
                  " (expected phi =", float(phi), ")"),
            if abs(pf_eigenvalue_numeric(M_fib) - float(phi)) > 1e-8
               then error("Fibonacci substitution PF eigenvalue != phi"))),

    /* ----------------------------------------------------------------
       Accumulate cross-section data for Section 20
       ---------------------------------------------------------------- */
    cross_section_data : endcons(
        [n_int,
         float(deficit_n * 180 / %pi),
         k_loop_close_n,
         float(Cn_gap),
         float(heat_trace_cycle(n_int, 1.0)),
         float(spectral_zeta_cycle(n_int, 2)),
         float(ngon_window_area(n_int, 1.0))],
        cross_section_data)

) $ /* end for N_GON */

/* ---------------------------------------------------------------------
    SECTION 10.  Q(sqrt(5))-arithmetic and units    [pentagon-specific,
                                                      run after sweep]
   --------------------------------------------------------------------- */

print("------------------------------------------------------------")$
print("[Pentagon-specific] Z[phi] integer arithmetic:")$

zphi_add(p, q) := [p[1] + q[1], p[2] + q[2]] $
zphi_neg(p)    := [-p[1], -p[2]] $
zphi_mul(p, q) := [p[1]*q[1] + p[2]*q[2],
                   p[1]*q[2] + p[2]*q[1] + p[2]*q[2]] $
zphi_to_sym(p) := p[1] + p[2]*phi $

e1 : [1, 0] $
e2 : [0, 1] $
phi_sq_pair : zphi_mul(e2, e2) $
if phi_sq_pair # [1, 1]
   then error("Z[phi] multiplication failed for phi^2")$
phi_cubed_pair : zphi_mul(phi_sq_pair, e2) $
if phi_cubed_pair # [1, 2]
   then error("Z[phi] multiplication failed for phi^3")$
print("  phi^2 = [1,1], phi^3 = [1,2] (Fibonacci pairs).")$

for nn:1 thru 7 do block(
    [p, expected],
    p : e2,
    for kk:2 thru nn do p : zphi_mul(p, e2),
    expected : [fib(nn-1), fib(nn)],
    if p # expected
       then error("Z[phi] form mismatch at n =", nn))$
print("  phi^n = [F_{n-1}, F_n] verified for n = 1..7.")$

norm_zphi(p) := p[1]^2 + p[1]*p[2] - p[2]^2 $
for nn:1 thru 5 do block(
    [p, n_int_l, n_sym],
    p : e2,
    for kk:2 thru nn do p : zphi_mul(p, e2),
    n_int_l : norm_zphi(p),
    n_sym : qnorm(zphi_to_sym(p)),
    if ratsimp(n_int_l - n_sym) # 0
       then error("Norm mismatch at phi^", nn))$
print("  Z[phi] norm formula N(a+b phi) = a^2 + ab - b^2 verified.")$

for nn:1 thru 7 do block(
    [p, n_int_l],
    p : e2,
    for kk:2 thru nn do p : zphi_mul(p, e2),
    n_int_l : norm_zphi(p),
    if n_int_l # (-1)^nn
       then error("Unit-norm test failed at n =", nn))$
print("  N(phi^n) = (-1)^n verified for n = 1..7 (phi^n a unit).")$
/* ---------------------------------------------------------------
    EXTENDED Z[phi] ARITHMETIC AND NUMBER THEORY
    --------------------------------------------------------------- */
print("------------------------------------------------------------")$
print("[Pentagon-specific] Extended Z[phi] number theory:")$
/* Division in Z[phi]: (a + b phi) / (c + d phi) = (a+b phi)(c+d psi) / N(c+d phi) */
zphi_conj(p) := [p[1] + p[2], -p[2]] $   /* conjugate: a+b phi -> a+b psi = (a+b) - b phi */
zphi_div(p, q) := block(
     [num, denom],
     num   : zphi_mul(p, zphi_conj(q)),
     denom : norm_zphi(q),
     if denom = 0 then error("Division by zero in Z[phi]")
     else [num[1]/denom, num[2]/denom]) $
/* GCD in Z[phi] via Euclidean algorithm (norm-based). */
zphi_gcd(p, q) := block(
     [a : p, b : q, r, na, nb, tmp],
     if norm_zphi(q) = 0 then return(p),
     if norm_zphi(p) = 0 then return(q),
     /* Ensure |N(a)| >= |N(b)| */
     if abs(norm_zphi(a)) < abs(norm_zphi(b)) then (tmp:a, a:b, b:tmp),
     for i : 1 thru 20 do block(
         [quot, rem],
         quot : zphi_div(a, b),
         quot : [round(quot[1]), round(quot[2])],
         rem  : zphi_add(a, zphi_neg(zphi_mul(b, quot))),
         if norm_zphi(rem) = 0 then return(b),
         a : b,
         b : rem),
     b) $
print("  Z[phi] conjugate: conj(phi) = [1,-1] (= 1 - phi = -psi):",
       zphi_conj(e2))$
if zphi_conj(e2) # [1, -1] then error("Z[phi] conjugate failed")$
print("  Z[phi] division: phi^2 / phi = phi:")$
block([q], q : zphi_div(phi_sq_pair, e2),
     print("    phi^2 / phi =", q, " (expected [0,1] = phi)"),
     if q # [0, 1] then error("Z[phi] division phi^2/phi failed"))$
print("  Z[phi] division: phi^3 / phi^2 = phi:")$
block([q], q : zphi_div(phi_cubed_pair, phi_sq_pair),
     print("    phi^3 / phi^2 =", q, " (expected [0,1] = phi)"),
     if q # [0, 1] then error("Z[phi] division phi^3/phi^2 failed"))$
/* Primes in Z[phi]: p splits if Legendre(5/p) = 1, inert if = -1, ramifies if p=5. */
print("  Prime splitting in Z[phi] = Z[(1+sqrt5)/2]:")$
for p_prime in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31] do block(
     [leg, split_type],
     leg : legendre5(p_prime),
     split_type : if leg = 0  then "ramified (p=5)"
                  else if leg = 1  then "splits (p = +-1 mod 5)"
                  else "inert   (p = +-2 mod 5)",
     print("    p =", p_prime, ": Legendre(5/p) =", leg, " =>", split_type))$
/* Units in Z[phi]: all phi^n for n in Z. */
print("  Units in Z[phi] (phi^n for n = -3..3):")$
for nn : -3 thru 3 do block(
     [u_sym, u_norm],
     u_sym  : ratsimp(phi^nn),
     u_norm : ratsimp(qnorm(u_sym)),
     print("    phi^", nn, " =", u_sym, " ~ ", float(u_sym),
           "  N =", u_norm))$
/* Continued fraction expansion of phi: [1; 1, 1, 1, ...] */
print("  Continued fraction convergents of phi:")$
block(
     [p_prev : 0, p_curr : 1, q_prev : 1, q_curr : 1],
     print("    h_0/k_0 = 1/1 =", float(1/1)),
     for i : 1 thru 10 do block(
         [p_next, q_next],
         p_next : p_curr + p_prev,
         q_next : q_curr + q_prev,
         print("    h_", i, "/k_", i, " =", p_next, "/", q_next,
               " =", float(p_next/q_next),
               "  error =", float(abs(p_next/q_next - phi))),
         p_prev : p_curr,  p_curr : p_next,
         q_prev : q_curr,  q_curr : q_next))$
/* Pisano period: Fibonacci sequence mod m. */
pisano_period(m) := block(
     [a : 0, b : 1, period : 0],
     for i : 1 thru 6*m do block(
         [c : mod(a + b, m)],
         a : b,  b : c,
         if a = 0 and b = 1 then (period : i, return(period))),
     period) $
print("  Pisano periods pi(m) for m = 2..15:")$
for m : 2 thru 15 do
     print("    pi(", m, ") =", pisano_period(m))$
/* Wall-Sun-Sun prime check (p^2 | F_{p - Legendre(5/p)}): */
print("  Wall-Sun-Sun prime check (F_{p - (5/p)} mod p^2) for small primes:")$
for p_prime in [5, 7, 11, 13, 17, 19, 23, 29] do block(
     [leg, idx, fib_val, wss],
     leg     : legendre5(p_prime),
     idx     : p_prime - leg,
     fib_val : fib(idx),
     wss     : mod(fib_val, p_prime^2),
     print("    p =", p_prime, ": F_{p-(5/p)} = F_", idx,
           " mod p^2 =", wss,
           if wss = 0 then " *** WALL-SUN-SUN PRIME ***" else ""))$

/* ---------------------------------------------------------------------
    SECTION 14.  Sweep summary table and final report
    --------------------------------------------------------------------- */

print("------------------------------------------------------------")$
print("SWEEP SUMMARY TABLE")$
print("  Range: n =", N_GON_MIN, "to", N_GON_MAX,
       "step", N_GON_STEP)$
print("------------------------------------------------------------")$
/* ---------------------------------------------------------------
    SECTION 20.  Cross-section correlation table
    --------------------------------------------------------------- */
print("------------------------------------------------------------")$
print("CROSS-SECTION CORRELATION TABLE")$
print("  Columns: n | deficit | k_close | C_n_gap | HeatTr(t=1) | SpZeta(s=2) | WinArea")$
print("------------------------------------------------------------")$
print(sconcat(
     pad_right("n",       4), " | ",
     pad_right("deficit", 9), " | ",
     pad_right("k_close", 8), " | ",
     pad_right("C_n_gap", 9), " | ",
     pad_right("HeatTr1", 9), " | ",
     pad_right("SpZeta2", 9), " | ",
     pad_right("WinArea", 9)))$
print(sconcat(
     "----+-----------+---------+----------+----------+----------+----------"))$
for row in cross_section_data do
     print(sconcat(
         pad_right(string(row[1]), 4), " | ",
         pad_right(string(float(row[2])), 9), " | ",
         pad_right(string(row[3]),  8), " | ",
         pad_right(string(float(row[4])), 9), " | ",
         pad_right(string(float(row[5])), 9), " | ",
         pad_right(string(float(row[6])), 9), " | ",
         pad_right(string(float(row[7])), 9)))$
print("------------------------------------------------------------")$
/* Correlation: deficit vs spectral gap */
print("Correlation analysis: deficit angle vs C_n spectral gap:")$
block(
     [deficits, gaps, n_pts, mean_d, mean_g, cov, var_d, var_g, corr],
     deficits : map(lambda([r], r[2]), cross_section_data),
     gaps     : map(lambda([r], r[4]), cross_section_data),
     n_pts    : length(deficits),
     mean_d   : float(lreduce("+", deficits, 0) / n_pts),
     mean_g   : float(lreduce("+", gaps,     0) / n_pts),
     cov      : float(lreduce("+",
                    map(lambda([d,g], (d - mean_d)*(g - mean_g)),
                        deficits, gaps), 0) / n_pts),
     var_d    : float(lreduce("+",
                    map(lambda([d], (d - mean_d)^2), deficits), 0) / n_pts),
     var_g    : float(lreduce("+",
                    map(lambda([g], (g - mean_g)^2), gaps),     0) / n_pts),
     corr     : if var_d * var_g > 0
                then cov / sqrt(var_d * var_g)
                else 0,
     print("  mean deficit =", mean_d, " deg"),
     print("  mean C_n gap =", mean_g),
     print("  Pearson r(deficit, C_n_gap) =", float(corr)))$
/* Correlation: k_close vs heat trace at t=1 */
print("Correlation: k_close vs heat trace Tr(e^{-L}) at t=1:")$
block(
     [kcs, hts, n_pts, mean_k, mean_h, cov, var_k, var_h, corr],
     kcs    : map(lambda([r], r[3]), cross_section_data),
     hts    : map(lambda([r], r[5]), cross_section_data),
     n_pts  : length(kcs),
     mean_k : float(lreduce("+", kcs, 0) / n_pts),
     mean_h : float(lreduce("+", hts, 0) / n_pts),
     cov    : float(lreduce("+",
                  map(lambda([k,h], (k - mean_k)*(h - mean_h)),
                      kcs, hts), 0) / n_pts),
     var_k  : float(lreduce("+",
                  map(lambda([k], (k - mean_k)^2), kcs), 0) / n_pts),
     var_h  : float(lreduce("+",
                  map(lambda([h], (h - mean_h)^2), hts), 0) / n_pts),
     corr   : if var_k * var_h > 0
              then cov / sqrt(var_k * var_h)
              else 0,
     print("  mean k_close =", mean_k),
     print("  mean HeatTr  =", mean_h),
     print("  Pearson r(k_close, HeatTr1) =", float(corr)))$

print(sconcat(
     pad_right("n",    4),  " | ",
     pad_right("theta(deg)", 12), " | ",
     pad_right("k_flat", 7), " | ",
     pad_right("deficit(deg)", 13), " | ",
     pad_right("k_close", 8), " | ",
     pad_right("turns", 6), " | ",
     pad_right("phi(n)", 7), " | ",
     pad_right("Area",   8), " | ",
     pad_right("R",      8), " | ",
     pad_right("OT_rules", 12), " | ",
     pad_right("C_n_gap", 9)))$
print(sconcat(
     "----+",
     "-------------+",
     "--------+",
     "--------------+",
     "---------+",
     "-------+",
     "--------+",
     "---------+",
     "---------+",
     "-------------+",
     "----------"))$
for row in ngon_summary_table do
     print(sconcat(
         pad_right(string(row[1]),  4),  " | ",
         pad_right(string(float(row[2])), 12), " | ",
         pad_right(string(row[3]),  7), " | ",
         pad_right(string(float(row[4])), 13), " | ",
         pad_right(string(row[5]),  8), " | ",
         pad_right(string(row[6]),  6), " | ",
         pad_right(string(row[7]),  7), " | ",
         pad_right(string(float(row[8])),  8), " | ",
         pad_right(string(float(row[9])),  8), " | ",
         pad_right(string(row[10]), 12), " | ",
         pad_right(string(float(row[11])), 9)))$
print("------------------------------------------------------------")$
print("Summary of verified identities and structural facts:")$
print(" [Generalized, swept over n =", N_GON_MIN, "to", N_GON_MAX, "]")$
print("  * Interior angle theta = (n-2) pi / n")$
print("  * Exterior angle = central angle = 2 pi / n")$
print("  * Flat-vertex count k_flat = floor(2 pi / theta)")$
print("  * Loop-closure index k_loop = 2n / gcd(2n, n-2)")$
print("  * R(2 pi/n)^n = I; det R = 1; tr R = 2 cos(2 pi / n)")$
print("  * R^k traces matched 2 cos(2 pi k/n) for k = 1..n")$
print("  * Primary basis vectors v[k] = (cos, sin)(pi/2 + 2 pi k/n)")$
print("  * Edge vectors e[k] computed and reported")$
print("  * Side length s = 2 sin(pi/n) (unit circumradius)")$
print("  * Z_n holonomy: order(g) = n/gcd(g,n)")$
print("  * Z_2 (spinor) and Z_n (anyonic) cases both verified")$
print("  * Alexander-Orbach: d_spec = 2 d_eff / d_w (n-independent)")$
print("  * charpoly(C_n) divides x^n - 1; C_n is orthogonal")$
print("  * Real-invariant 2D subspaces V_k of C_n enumerated")$
print("  * Projectors P_diag, P_perp idempotent, orthogonal, complete")$
print("  * dim(E_perp) = n - 2 (after taking V_1 as physical plane)")$
print("  * Outer-totalistic n-neighbor binary rule count = 2^(2(n+1))")$
print("  * Diagonals d[k] = sin(k pi/n)/sin(pi/n)")$
print("  * Area A = (n/4) cot(pi/n)")$
print("  * Dihedral D_n: |D_n| = 2n, S R S = R^{-1}")$
print("  * D_n irrep dim^2 sum = 2n (case n even / n odd)")$
print("  * Star K_{1,n} Laplacian spectrum = {0, 1 x n, n+1}")$
print("  * Cycle C_n Laplacian eigenvalues = 2 - 2 cos(2 pi k/n)")$
print(" [Pentagon-specific benchmarks verified]")$
print("  * phi^2 = phi + 1; phi*(phi-1) = 1; 1/phi = phi - 1")$
print("  * Fibonacci recursion phi^n = F_n phi + F_{n-1}")$
print("  * Binet formula F_n = (phi^n - psi^n)/sqrt(5)")$
print("  * Pentagon diagonal/side = phi")$
print("  * Pentagram self-similarity: 1/phi^2 = 2 - phi")$
print("  * Z[phi] arithmetic, N(phi^n) = (-1)^n")$
print(" [Extended sections verified]")$
print("  * Alexander-Orbach phase diagram (d_eff x d_w grid)")$
print("  * Anomalous diffusion classification by d_w")$
print("  * Crossover scale t_cross from UV to IR regime")$
print("  * Multifractal tau(q) spectrum (monofractal baseline)")$
print("  * Flat U(1) holonomy phases and Aharonov-Bohm phases")$
print("  * Monodromy matrices and loop orders for n=3,4,5,6")$
print("  * Berry phases and Chern numbers for Z_n sectors")$
print("  * Inflation/substitution matrices and PF eigenvalues")$
print("  * Pentagon inflation ratio = phi^2 confirmed")$
print("  * Dedekind zeta and Dirichlet L(s,chi_5) partial sums")$
print("  * Spectral zeta of C_n at s=1,2,3")$
print("  * Heat trace Tr(e^{-tL}) for cycle and star graphs")$
print("  * Effective spectral dimension from heat-trace slope")$
print("  * Acceptance-window geometry and density")$
print("  * Pentagon window inflation: radius -> 1/phi")$
print("  * Fibonacci word: length = F_k, a/b ratio -> phi")$
print("  * Fibonacci substitution charpoly = x^2 - x - 1")$
print("  * Thue-Morse and Rudin-Shapiro sequences generated")$
print("  * Z[phi] division, conjugate, GCD")$
print("  * Prime splitting in Z[phi] by Legendre symbol")$
print("  * Continued fraction convergents of phi")$
print("  * Pisano periods pi(m) for m=2..15")$
print("  * Wall-Sun-Sun prime check for small primes")$
print("  * Cross-section correlation table (20 columns)")$
print("  * Pearson correlations: deficit/gap, k_close/HeatTr")$
print("============================================================")$
print("analysis.mac : all checks passed for n =",
       N_GON_MIN, "to", N_GON_MAX, "step", N_GON_STEP)$
print("============================================================")$