Scale up the empirical figures using the fast factorization

.github/workflows/ci.yml

@@ -13,6 +13,6 @@ jobs:
         with:
           python-version: "3.12"
       - name: Install test dependencies
-        run: pip install pytest
+        run: pip install -r requirements.txt pytest
       - name: Run tests
         run: pytest tests/ -v

README.md

@@ -121,16 +121,34 @@ python propagation.py
 ### Distribution of `gcd(p-1, q-1)` (`gcd_distribution_theory.py`)
 
 Empirical verification of Theorem C. Computes the divisibility rate
-`Pr(l | gcd(p-1, q-1))` over all distinct odd-prime pairs up to a
-configurable cutoff and compares it to the Dirichlet prediction
-`1/(l-1)^2`. Also prints the empirical mean gcd alongside the
-asymptotic estimate `A log X` with `A = 315 zeta(3) / (2 pi^4) ~ 1.94`.
+`Pr(l | gcd(p-1, q-1))` over distinct odd-prime pairs up to a
+configurable cutoff (`10^5` by default, switching to a fixed-seed
+random sample of 2M pairs above the exhaustive threshold) and compares
+it to the Dirichlet prediction `1/(l-1)^2`. At the `10^5` cutoff the
+relative errors are around a percent. Also prints the empirical mean
+gcd alongside the asymptotic estimate `A log X` with
+`A = 315 zeta(3) / (2 pi^4) ~ 1.94`.
 
 ```bash
 python gcd_distribution_theory.py
 # -> gcd_distribution.png
 ```
 
+### Collapse at scale (`collapse_at_scale.py`)
+
+The full-range scan of `C(n)` to `n = 10^6`, using a smallest-prime-factor
+sieve (one pass factors every n at once, so the scan takes seconds).
+Produces a 2-panel figure: the density of `(n, C)` over all composites
+with the running-maximum "collapse records" overlaid and labeled, and a
+quantile fan (median / 90th / 99th / max of `C` by geometric window)
+showing that collapse is a tail phenomenon. Prints the by-kind summary
+table and the top records with factorizations.
+
+```bash
+python collapse_at_scale.py
+# -> collapse_at_scale.png
+```
+
 ### PDF of the theorem note (`render_theorem.ps1`)
 
 ```powershell
@@ -159,6 +177,21 @@ products of small primes whose `(p_i - 1)` shares many common factors
 The Hardy-Ramanujan number `1729 = 7 * 13 * 19` shows up near the top:
 it is also a Carmichael number, and large `C` is part of why.
 
+Scaling up with `collapse_at_scale.py` over `n in [2, 10^6]`:
+
+| kind         | count  | mean C | median C | max C | fraction with C=1 |
+| ------------ | ------ | ------ | -------- | ----- | ----------------- |
+| prime        | 78498  | 1.00   | 1        | 1     | 1.000             |
+| prime_power  | 236    | 1.07   | 1        | 2     | 0.928             |
+| composite    | 921265 | 36.22  | 8        | 10368 | 0.045             |
+
+The record holder below `10^6` is `959595 = 3 * 5 * 7 * 13 * 19 * 37`
+with `C = 10368` -- the same story as `1365`, two primes deeper: the
+totients `2, 4, 6, 12, 18, 36` all divide each other's lattice. The
+median composite has `C = 8` while the maximum is `10368`, i.e. the
+mean is dragged by a thin tail of heavily-shared-structure numbers;
+the quantile fan in `collapse_at_scale.png` makes this visible.
+
 ## Cryptographic interpretation
 
 For `n = p*q` (an RSA-shaped modulus), the order of an arbitrary unit
@@ -200,6 +233,7 @@ identity at runtime, so they double as tests for the theorem in
 | `wedge_envelopes.py`       | Algebraic envelope visualization                        |
 | `propagation.py`           | Iterative demo of the collapse propagation theorem      |
 | `gcd_distribution_theory.py` | Empirical verification of Theorem C                   |
+| `collapse_at_scale.py`     | Sieve-based scan of `C(n)` to `10^6`, records + quantiles |
 | `THEOREM.md`               | Wedge, propagation, and Dirichlet-density identities    |
 | `render_theorem.ps1`       | Pandoc helper that renders `THEOREM.md` to `theorem.pdf`|
 | `requirements.txt`         | `matplotlib` (pulls in numpy)                           |

THEOREM.md

@@ -287,6 +287,7 @@ multiplicative law shown at three scales:
 | `order_distributions.png`  | The parallel-cycle picture inside specific $n$, with explicit invariant factors. |
 | `propagation.png`          | Theorem 2 made visible: stacked $\log_2 \gcd$ contributions per prime, plus the empirical density of $\gcd(p - 1, q - 1)$. |
 | `gcd_distribution.png`     | Theorem C made visible: empirical vs Dirichlet-predicted divisibility rates, plus the $1365$ tie-back. |
+| `collapse_at_scale.png`    | The whole story to $n = 10^6$: collapse records (the $1365$-style extremes, deeper) and the quantile fan showing collapse is a tail phenomenon. |
 
 ## 10. Honest scope
 

collapse_at_scale.py

@@ -0,0 +1,211 @@
+#!/usr/bin/env python3
+"""
+collapse_at_scale.py
+
+The collapse index C(n) = phi(n)/lambda(n) scanned to n = 10^6.
+
+The per-n factorization in lambda_ratio_explorer.py is the right tool for
+individual (possibly huge) n, but a full-range scan is faster with a
+smallest-prime-factor sieve: one pass builds the factorization of every
+n <= N, from which phi and lambda follow directly.
+
+Generates a 2-panel figure: collapse_at_scale.png
+  1. Density of (n, C(n)) over all composites on log-log axes, with the
+     running-maximum "collapse records" overlaid and labeled. The record
+     holders are exactly the products of small primes with heavily shared
+     totient structure that Theorem C predicts.
+  2. Quantile fan of C(n) by scale: median, 90th, 99th percentile, and max
+     of C over geometric windows in n. The median stays in single digits
+     while the extremes explode -- collapse is a tail phenomenon.
+
+Also prints the by-kind summary table (the 10^6 version of the README
+table) and the top collapse records with factorizations.
+"""
+
+from __future__ import annotations
+
+import math
+from array import array
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from lambda_ratio_explorer import factorize
+
+N_MAX = 1_000_000
+N_QUANTILE_BINS = 40
+
+# Repo palette: blues carry magnitude (sequential, light -> dark),
+# orange marks the record trace, grays are ink.
+BLUE_SEQ = ["#74c0fc", "#339af0", "#1971c2", "#0b4a8c"]
+RECORD_COLOR = "#e8590c"
+INK = "#495057"
+
+
+def spf_sieve(limit: int) -> array:
+    """Smallest prime factor of every n <= limit."""
+    spf = array("l", range(limit + 1))
+    for i in range(2, int(limit ** 0.5) + 1):
+        if spf[i] == i:
+            for j in range(i * i, limit + 1, i):
+                if spf[j] == j:
+                    spf[j] = i
+    return spf
+
+
+def scan(limit: int) -> tuple[np.ndarray, np.ndarray]:
+    """Return (C, kind_code) arrays indexed by n.
+
+    kind_code: 0 = unit/unused, 1 = prime, 2 = prime power, 3 = composite.
+    """
+    spf = spf_sieve(limit)
+    C = np.zeros(limit + 1, dtype=np.int64)
+    kind_code = np.zeros(limit + 1, dtype=np.int8)
+    C[1] = 1
+
+    for n in range(2, limit + 1):
+        m = n
+        phi = 1
+        lam = 1
+        distinct = 0
+        while m > 1:
+            p = spf[m]
+            k = 0
+            while m % p == 0:
+                m //= p
+                k += 1
+            distinct += 1
+            phi *= (p - 1) * p ** (k - 1)
+            if p == 2 and k >= 3:
+                lam_pk = 2 ** (k - 2)
+            else:
+                lam_pk = (p - 1) * p ** (k - 1)
+            lam = lam * lam_pk // math.gcd(lam, lam_pk)
+            if distinct == 1 and m == 1:
+                kind_code[n] = 1 if k == 1 else 2
+        if kind_code[n] == 0:
+            kind_code[n] = 3
+        C[n] = phi // lam
+    return C, kind_code
+
+
+def collapse_records(C: np.ndarray) -> list[tuple[int, int]]:
+    """(n, C(n)) points where C sets a new running maximum."""
+    records = []
+    best = 0
+    for n in range(2, len(C)):
+        if C[n] > best:
+            best = int(C[n])
+            records.append((n, best))
+    return records
+
+
+def print_summary(C: np.ndarray, kind_code: np.ndarray) -> None:
+    names = {1: "prime", 2: "prime_power", 3: "composite"}
+    print(f"\nCollapse index C(n) = phi(n)/lambda(n) over n in [2, {len(C) - 1}]\n")
+    print(f"  {'kind':<12} {'count':>8} {'mean C':>9} {'median C':>9} {'max C':>8} {'frac C=1':>9}")
+    for code, name in names.items():
+        vals = C[2:][kind_code[2:] == code]
+        print(f"  {name:<12} {len(vals):>8} {vals.mean():>9.2f} {int(np.median(vals)):>9} "
+              f"{vals.max():>8} {(vals == 1).mean():>9.3f}")
+
+    print("\nTop collapse records:")
+    records = collapse_records(C)
+    for n, c in records[-10:]:
+        factors = " * ".join(
+            (f"{p}^{k}" if k > 1 else str(p)) for p, k in factorize(n).items()
+        )
+        print(f"  C({n:>7}) = {c:>5}   {n} = {factors}")
+
+
+def panel_density(ax: plt.Axes, C: np.ndarray, kind_code: np.ndarray) -> None:
+    comp = np.flatnonzero(kind_code == 3)
+    x = np.log10(comp.astype(np.float64))
+    y = np.log2(C[comp].astype(np.float64))
+
+    hb = ax.hexbin(x, y, gridsize=60, cmap="Blues", bins="log", mincnt=1,
+                   linewidths=0.1)
+    cb = plt.colorbar(hb, ax=ax, pad=0.01)
+    cb.set_label("composites per cell (log scale)", fontsize=8.5)
+
+    records = collapse_records(C)
+    rx = [math.log10(n) for n, _ in records]
+    ry = [math.log2(c) for _, c in records]
+    ax.plot(rx, ry, color=RECORD_COLOR, lw=1.4, zorder=3,
+            drawstyle="steps-post", label="running max of $C$")
+    ax.scatter(rx, ry, color=RECORD_COLOR, s=22, zorder=4,
+               edgecolor="white", linewidth=0.6)
+    # Label a spread of records: the final one plus two mid-scale ones,
+    # so annotations never pile up in the top-right corner.
+    labeled = {records[-1]}
+    for frac in (0.5, 0.78):
+        target = frac * math.log10(records[-1][0])
+        labeled.add(min(records, key=lambda r: abs(math.log10(r[0]) - target)))
+    for n, c in labeled:
+        ax.annotate(f"$C({n}) = {c}$", (math.log10(n), math.log2(c)),
+                    xytext=(-8, 6), textcoords="offset points",
+                    ha="right", fontsize=8.5, color=INK)
+
+    xticks = range(1, 7)
+    ax.set_xticks(list(xticks))
+    ax.set_xticklabels([f"$10^{k}$" for k in xticks])
+    yticks = range(0, int(max(ry)) + 2, 2)
+    ax.set_yticks(list(yticks))
+    ax.set_yticklabels([f"{2 ** k}" for k in yticks])
+    ax.set_xlabel("$n$")
+    ax.set_ylabel("$C(n)$")
+    ax.set_title("Collapse of composites to $n = 10^6$\n"
+                 "density of $(n, C)$ with record collapses overlaid")
+    ax.grid(True, alpha=0.15)
+    ax.legend(loc="upper left", fontsize="small")
+
+
+def panel_quantile_fan(ax: plt.Axes, C: np.ndarray, kind_code: np.ndarray) -> None:
+    comp = np.flatnonzero(kind_code == 3)
+    Cc = C[comp].astype(np.float64)
+    edges = np.geomspace(comp[0], len(C) - 1, N_QUANTILE_BINS + 1)
+
+    quantiles = [(0.50, "median", BLUE_SEQ[0]),
+                 (0.90, "90th pct", BLUE_SEQ[1]),
+                 (0.99, "99th pct", BLUE_SEQ[2]),
+                 (1.00, "max", BLUE_SEQ[3])]
+    centers = np.sqrt(edges[:-1] * edges[1:])
+    for q, label, color in quantiles:
+        ys = []
+        for lo, hi in zip(edges[:-1], edges[1:]):
+            window = Cc[(comp >= lo) & (comp < hi)]
+            ys.append(np.quantile(window, q) if len(window) else np.nan)
+        ax.plot(centers, ys, color=color, lw=1.8)
+        ax.annotate(label, (centers[-1], ys[-1]), xytext=(6, 0),
+                    textcoords="offset points", va="center",
+                    fontsize=8.5, color=INK)
+
+    ax.set_xscale("log")
+    ax.set_yscale("log", base=2)
+    ax.set_xlabel("$n$ (geometric windows)")
+    ax.set_ylabel("$C(n)$ quantile within window")
+    ax.set_title("Collapse is a tail phenomenon\n"
+                 "the median composite barely collapses; the extremes explode")
+    ax.grid(True, which="both", alpha=0.15)
+    ax.set_xlim(centers[0], centers[-1] * 3.2)
+
+
+def main() -> None:
+    print(f"sieving and scanning to {N_MAX} ...")
+    C, kind_code = scan(N_MAX)
+    print_summary(C, kind_code)
+
+    fig, axes = plt.subplots(1, 2, figsize=(15, 6.2))
+    fig.suptitle("Collapse index $C(n) = \\varphi(n)/\\lambda(n)$ at scale",
+                 fontsize=13, fontweight="bold")
+    panel_density(axes[0], C, kind_code)
+    panel_quantile_fan(axes[1], C, kind_code)
+
+    fig.tight_layout()
+    out = "collapse_at_scale.png"
+    fig.savefig(out, dpi=160)
+    print(f"\nWrote {out}")
+
+
+if __name__ == "__main__":
+    main()

gcd_distribution_theory.py

@@ -12,8 +12,10 @@
 This script:
   1. Sieves odd primes up to a configurable limit.
   2. For l in {2, 3, 5, 7, 11, 13, 17, 19, 23} computes empirical
-     Pr( l | gcd(p - 1, q - 1) ) over all distinct odd prime pairs and
-     compares to the Dirichlet prediction 1/(l-1)^2.
+     Pr( l | gcd(p - 1, q - 1) ) over distinct odd prime pairs and
+     compares to the Dirichlet prediction 1/(l-1)^2. Below a size
+     threshold all pairs are enumerated; above it a fixed-seed random
+     sample of pairs is used so the cutoff can be pushed to 10^5.
   3. Computes empirical E[gcd(p - 1, q - 1)] and the truncated heuristic
      sum_{d <= D} 1/phi(d), which grows like A log D with
      A = 315 zeta(3) / (2 pi^4).
@@ -30,6 +32,7 @@
 from __future__ import annotations
 
 import math
+import random
 from typing import Iterable
 
 import matplotlib.pyplot as plt
@@ -41,9 +44,15 @@
 )
 
 
-PRIME_LIMIT = 2000
+PRIME_LIMIT = 100_000
 TARGET_PRIMES: list[int] = [2, 3, 5, 7, 11, 13, 17, 19, 23]
 
+# Above this many pairs, sample instead of enumerating. The seed is fixed
+# so the table and figure are reproducible run to run.
+MAX_EXHAUSTIVE_PAIRS = 2_000_000
+SAMPLE_PAIRS = 2_000_000
+SAMPLE_SEED = 20260711
+
 # Mean-totient-reciprocal constant A = 315 zeta(3) / (2 pi^4),
 # governing the asymptotic sum_{d <= X} 1/phi(d) ~ A log X.
 ASYMPTOTIC_CONSTANT_A = 1.9436
@@ -61,21 +70,37 @@ def primes_up_to(limit: int) -> list[int]:
     return [i for i, b in enumerate(sieve) if b]
 
 
+def _pair_stream(primes: list[int]) -> Iterable[tuple[int, int]]:
+    """Distinct prime pairs: exhaustive when small, sampled when large."""
+    n = len(primes)
+    if n * (n - 1) // 2 <= MAX_EXHAUSTIVE_PAIRS:
+        for i, p in enumerate(primes):
+            for q in primes[i + 1:]:
+                yield p, q
+        return
+    rng = random.Random(SAMPLE_SEED)
+    for _ in range(SAMPLE_PAIRS):
+        i = rng.randrange(n)
+        j = rng.randrange(n - 1)
+        if j >= i:
+            j += 1
+        yield primes[i], primes[j]
+
+
 def empirical_divisibility_rates(
     primes: list[int], target_ells: Iterable[int]
 ) -> tuple[dict[int, float], int, float]:
-    """Return (empirical Pr(l | gcd) per l, total pair count, empirical mean gcd)."""
+    """Return (empirical Pr(l | gcd) per l, pair count used, empirical mean gcd)."""
     counts = {l: 0 for l in target_ells}
     pair_total = 0
     gcd_sum = 0
-    for i, p in enumerate(primes):
-        for q in primes[i + 1:]:
-            g = math.gcd(p - 1, q - 1)
-            pair_total += 1
-            gcd_sum += g
-            for l in counts:
-                if g % l == 0:
-                    counts[l] += 1
+    for p, q in _pair_stream(primes):
+        g = math.gcd(p - 1, q - 1)
+        pair_total += 1
+        gcd_sum += g
+        for l in counts:
+            if g % l == 0:
+                counts[l] += 1
     rates = {l: counts[l] / pair_total for l in counts}
     mean_gcd = gcd_sum / pair_total
     return rates, pair_total, mean_gcd