Frequent strong liars among witness numbers (with top 100 results up to 115,471)

liars.csv

liars.py

import csv
import logging
import time

from argparse import ArgumentParser
from typing import List, Optional, Tuple
from collections import Counter


class Suspect:
    """ An integer accused of being prime. """
    __slots__ = ["n", "s", "d"]
    def __init__(self, n: int):
        s, d = prepare_suspect(n)
        self.n = n
        self.s = s
        self.d = d


def prepare_suspect(n: int) -> Tuple[int, int]:
    """ Find integers s and d such that n = 2^s * d + 1 and d is odd. """
    s = 1
    d, r = divmod(n, 2)
    if not r:
        raise ValueError("`n` must be odd.")

    while True:
        q, r = divmod(d, 2)
        if r:
            break
        else:
            s += 1
            d = q
        
    return s, d


def is_prime(s: Suspect, witnesses: Optional[List[int]] = None) -> bool:
    """ Check if a suspect is prime based on a set of witnesses. """
    if not witnesses:
        witnesses = call_star_witnesses(s.n)
    return all(witness_statement(w, s) == True for w in witnesses)


def witness_statement(witness: int, suspect: Suspect) -> bool:
    """ Do any Miller-Rabin congruence relations hold? """
    a = witness
    n = suspect.n
    s = suspect.s
    d = suspect.d
    x = pow(a, d, n)
    if x in {1, n - 1}:
        return True
    while (s := s - 1):
        x = pow(x, 2, n)
        if x == (n - 1):
            return True
    return False


MAX_SUSPECT = 1_122_004_669_633

def call_star_witnesses(n: int):
    """
    Find sets of witnesses whose agreement is sufficient to declare primality.
    Presented by Matt Parker on Numberphile: https://youtu.be/_MscGSN5J6o
    """
    if n < 1_373_653:
        return [2, 3]
    elif n < 25_326_001:
        return [2, 3, 5]
    elif n < MAX_SUSPECT:
        return [2, 13, 23, 1_662_803]
    else:
        raise ValueError(f"Don't know star witnesses for x >= {MAX_SUSPECT}.")


if __name__ == "__main__":
    parser = ArgumentParser()
    parser.add_argument("n", type=int, nargs="?", default=100)
    parser.add_argument("--top", type=int, default=5)
    parser.add_argument("--verbose", action="store_true", default=False)
    parser.add_argument("--output", type=str, default="")
    args = parser.parse_args()

    logging.basicConfig(level=logging.DEBUG if args.verbose else logging.INFO)
    
    start = lap_start = time.perf_counter()
    
    lies: Counter[int] = Counter()
    try:
        # Test all "suspects" up to given `n`
        for n in range(5, args.n, 2):
            s = Suspect(n)

            # Get true verdict from star witnesses
            star_witnesses = call_star_witnesses(n)
            truth = is_prime(s, witnesses=star_witnesses)
            if truth:
                logging.debug(f"{n} is prime")
            else:
                logging.debug(f"{n} is not prime")
            
            # Get statements from all potential witnesses and count lies
            potential_witnesses = range(2, n)
            for w in potential_witnesses:
                if witness_statement(w, s) != truth:
                    lies[w] += 1
            
            # Output some progress info every now and again
            if (n + 1) % 1000 == 0:
                lap_stop = time.perf_counter()
                logging.info(f"Done {n + 1}/{args.n} ... ({lap_stop - lap_start:.2f}s)")
                lap_start = lap_stop
            
    except KeyboardInterrupt:
        pass # Finish reporting statistics before exiting

    stop = time.perf_counter()
    
    logging.info(f"Done ({stop - start:.2f}s total)")
    logging.info(f"The top {args.top} liars up to {n} are:")
    logging.info(lies.most_common(args.top))

    if args.output:
        with open(args.output, "w", newline="") as f:
            w = csv.DictWriter(f, fieldnames=["witness", "lies"])
            w.writeheader()
            for witness, lies in lies.most_common():
                w.writerow({ "witness" : witness, "lies" : lies })
            logging.info(f"Full results written to '{args.output}'.")