diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 64a6db1..4be0aed 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/README.md b/README.md index f8c06e8..164ccf9 100644 --- a/README.md +++ b/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) | diff --git a/THEOREM.md b/THEOREM.md index 2e62a2d..5775d6a 100644 --- a/THEOREM.md +++ b/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 diff --git a/collapse_at_scale.png b/collapse_at_scale.png new file mode 100644 index 0000000..7a08ad2 Binary files /dev/null and b/collapse_at_scale.png differ diff --git a/collapse_at_scale.py b/collapse_at_scale.py new file mode 100644 index 0000000..b8fecf4 --- /dev/null +++ b/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() diff --git a/gcd_distribution.png b/gcd_distribution.png index 60383c1..4b3d583 100644 Binary files a/gcd_distribution.png and b/gcd_distribution.png differ diff --git a/gcd_distribution_theory.py b/gcd_distribution_theory.py index 3b4f2b5..6128ec0 100644 --- a/gcd_distribution_theory.py +++ b/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 diff --git a/tests/test_core.py b/tests/test_core.py index 50b0d20..7ce67c6 100644 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -206,6 +206,18 @@ def test_matches_lambda_criterion(self, n): assert is_carmichael(n) == expected +class TestCollapseAtScaleScan: + def test_sieve_scan_matches_library(self): + pytest.importorskip("matplotlib") + from collapse_at_scale import scan + + C, kind_code = scan(2000) + kind_names = {1: "prime", 2: "prime_power", 3: "composite"} + for n in range(2, 2001): + assert C[n] == collapse_index(n) + assert kind_names[kind_code[n]] == kind(n) + + class TestKind: def test_classification(self): assert kind(7) == "prime"