"""Reproduce the elementary numerical claims in the primality-testing post."""

from __future__ import annotations

import math


def is_prime(n: int) -> bool:
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    return all(n % divisor for divisor in range(3, math.isqrt(n) + 1, 2))


def base_two_fermat_pseudoprimes(limit: int) -> list[int]:
    return [
        n
        for n in range(3, limit, 2)
        if not is_prime(n) and math.gcd(2, n) == 1 and pow(2, n - 1, n) == 1
    ]


def main() -> None:
    pseudoprimes = base_two_fermat_pseudoprimes(1_000_000)

    print(f"2^560 mod 561 = {pow(2, 560, 561)}")
    print(f"5^560 mod 561 = {pow(5, 560, 561)}")
    print(f"sqrt(10^10) = {math.isqrt(10**10):,}")
    print("sqrt(10^300) = 10^150")
    print("10^150 divisions / 10^12 per second = 10^138 seconds")
    print(f"base-2 pseudoprimes below 1,000: {[n for n in pseudoprimes if n < 1_000]}")
    print(f"base-2 pseudoprimes below 1,000,000: {len(pseudoprimes)}")
    print(f"341 chain: 2^85 mod 341 = {pow(2, 85, 341)}")
    print(f"then 32^2 mod 341 = {pow(32, 2, 341)}")
    print(f"4^-40 = {4**-40:.16g} = 2^-80")
    expected_odd_candidates = math.log(2**1024) / 2
    print(f"ln(2^1024) / 2 = {expected_odd_candidates:.6f}")

    assert 561 == 3 * 11 * 17
    assert all(560 % divisor == 0 for divisor in (2, 10, 16))
    assert [n for n in pseudoprimes if n < 1_000] == [341, 561, 645]
    assert len(pseudoprimes) == 245


if __name__ == "__main__":
    main()