Pratt Clifford Zeta Golomb spigot

cform_order_twin_spigot.py

#!/usr/bin/env python3
from __future__ import annotations

import argparse
import json
from functools import lru_cache
from typing import Dict, List, Tuple

import mpmath as mp
import sympy as sp


def smooth_products_with_factors(bases: List[int], bound: int) -> Dict[int, Dict[int, int]]:
    out: Dict[int, Dict[int, int]] = {1: {}}
    bs = sorted(set(bases))
    for p in bs:
        cur_items = list(out.items())
        for m, fac in cur_items:
            v = m
            e = 0
            while v <= bound // p:
                v *= p
                e += 1
                if v not in out:
                    nf = dict(fac)
                    nf[p] = e
                    out[v] = nf
    return out


def discover_primes_by_cform_closure(pmax: int) -> Tuple[Dict[int, Dict[int, int]], int]:
    # Closure rule only: p is admitted iff p = 1 + M, where M is smooth over known primes.
    info: Dict[int, Dict[int, int]] = {2: {}}
    rounds = 0
    while True:
        rounds += 1
        known = sorted(info.keys())
        smooth = smooth_products_with_factors(known, pmax - 1)
        new = 0
        for m, fac in smooth.items():
            cand = m + 1
            if cand > pmax or cand in info:
                continue
            if sp.isprime(cand):
                info[cand] = dict(fac)
                new += 1
        if new == 0:
            break
    return info, rounds


def cform_ordered_primes(pmax: int) -> Tuple[List[int], Dict[int, Dict[int, int]], int]:
    info, rounds = discover_primes_by_cform_closure(pmax)

    @lru_cache(None)
    def cform_radical(p: int) -> str:
        if p == 2:
            return "()"
        kids = sorted(cform_radical(q) for q in sorted(info[p].keys()))
        return "(" + "".join(kids) + ")"

    ps = sorted(info.keys(), key=lambda p: (len(cform_radical(p)), cform_radical(p), p))
    return ps, info, rounds


def mul_inv_geom(poly: List[int], k: int, a: int, deg: int) -> List[int]:
    # Multiply by 1/(1-a x^k) truncated to degree deg.
    out = poly[:]
    for d in range(k, deg + 1):
        out[d] += a * out[d - k]
    return out


def mul_one_minus(poly: List[int], k: int, a: int, deg: int) -> List[int]:
    # Multiply by (1-a x^k) truncated to degree deg.
    out = poly[:]
    for d in range(deg, k - 1, -1):
        out[d] -= a * out[d - k]
    return out


def build_euler_series(primes: List[int], deg: int) -> List[int]:
    f = [0] * (deg + 1)
    f[0] = 1
    for i, p in enumerate(primes, start=1):
        if i > deg:
            break
        f = mul_inv_geom(f, i, p, deg)
    return f


def recover_prefix_from_series(series: List[int], ordered_primes: List[int], m: int) -> List[int]:
    deg = len(series) - 1
    rec = []
    for n in range(1, m + 1):
        g = series[:]
        for i in range(1, n):
            g = mul_one_minus(g, i, ordered_primes[i - 1], deg)
        rec.append(g[n])
    return rec


def numeric_limit_estimate(ordered_primes: List[int], n: int, n_terms: int, x: mp.mpf) -> mp.mpf:
    # U_n(x)=prod_{k=n..n_terms}(1-p_k x^k)^-1, then p_n ~ (U_n-1)/x^n as x->0
    u = mp.mpf("1.0")
    for k in range(n, n_terms + 1):
        p = ordered_primes[k - 1]
        u *= 1 / (1 - p * (x ** k))
    return (u - 1) / (x ** n)


def find_twins(values: List[int]) -> Tuple[List[Tuple[int, int]], List[Tuple[int, int, int]]]:
    s = set(values)
    by_value = sorted((p, p + 2) for p in s if (p + 2) in s)
    by_order = []
    for i in range(len(values) - 1):
        if values[i + 1] - values[i] == 2:
            by_order.append((i + 1, values[i], values[i + 1]))
    return by_value, by_order


def run(pmax: int, m: int, x_str: str) -> dict:
    ordered, info, rounds = cform_ordered_primes(pmax)
    if len(ordered) < m:
        raise ValueError(f"need at least {m} primes in c-form order, only {len(ordered)} up to pmax={pmax}")

    @lru_cache(None)
    def cform_radical(p: int) -> str:
        if p == 2:
            return "()"
        kids = sorted(cform_radical(q) for q in sorted(info[p].keys()))
        return "(" + "".join(kids) + ")"

    deg = m + 8
    series = build_euler_series(ordered, deg)
    recovered = recover_prefix_from_series(series, ordered, m)

    ok = recovered == ordered[:m]

    # Numerical limit sanity checks for first few terms.
    mp.mp.dps = 100
    x = mp.mpf(x_str)
    numeric = []
    for n in range(1, min(8, m) + 1):
        est = numeric_limit_estimate(ordered, n, min(len(ordered), m + 12), x)
        numeric.append({
            "n": n,
            "target": ordered[n - 1],
            "estimate": str(est),
            "rounded": int(mp.nint(est)),
        })

    twins_value, twins_order = find_twins(recovered)

    growth = []
    for mm in [20, 40, 60, 80, 100]:
        if mm > m:
            continue
        tval, _ = find_twins(recovered[:mm])
        growth.append({"m": mm, "twin_pairs_in_prefix": len(tval)})

    head = [{"n": i + 1, "p": p, "shape": cform_radical(p)} for i, p in enumerate(ordered[:min(m, 30)])]
    return {
        "pmax": pmax,
        "m": m,
        "closure_rounds": rounds,
        "n_primes_discovered": len(ordered),
        "order_key": "(len(radical_cform), radical_cform_lex, prime)",
        "recovery_exact": ok,
        "ordered_head": head,
        "recovered_head": recovered[:min(m, 30)],
        "numeric_limit_checks": numeric,
        "twins_by_value_in_prefix": twins_value,
        "twins_adjacent_in_cform_order": twins_order,
        "twin_growth": growth,
    }


def main() -> None:
    ap = argparse.ArgumentParser(description="Twin-prime spigot via c-form ordered Euler product")
    ap.add_argument("--pmax", type=int, default=5000)
    ap.add_argument("--m", type=int, default=120)
    ap.add_argument("--x", type=str, default="1e-4")
    ap.add_argument("--out", type=str, default="clifford_zeta_lab/cform_twin_spigot_results.json")
    args = ap.parse_args()

    out = run(args.pmax, args.m, args.x)
    with open(args.out, "w", encoding="ascii") as f:
        json.dump(out, f, indent=2)

    print(json.dumps({
        "out": args.out,
        "recovery_exact": out["recovery_exact"],
        "n_twins_by_value": len(out["twins_by_value_in_prefix"]),
        "n_twins_adjacent_cform_order": len(out["twins_adjacent_in_cform_order"]),
    }, indent=2))


if __name__ == "__main__":
    main()