#!/usr/bin/env python3
import random
import time
import math
import sys

def rabin_miller_deterministic(n, cores=1, worker=0):
    def test_once(a,d,r,n):
        x = pow(a, d, n)
        if x == 1 or x == n-1: return True
        for i in range(1,r):
            x = (x*x) % n
            if x == n-1: return True
        return False

    # Take care of special cases
    if n == 2 or n == 3: return True
    if n == 1 or n % 2 == 0: return False

    # Find r, d such that n = 2**r * d + 1, with d odd
    d = n - 1
    r = 0
    while d % 2 == 0:
        d = d // 2
        r += 1

    ul = min(n-2, math.floor(2 * math.log(n)**2))
    if worker == 0: print("ul =", ul)

    for a in range(2, ul+1):
        if cores == 1 or a % cores == worker:
            if test_once(a,d,r,n) == False:
                return False
    return True

def rabin_miller(n, certainty=128):
    """Determine if a positive integer is (probably) a prime.

    Positional argument:

        n   --      The positive integer to check for primality.

    Keyword argument:

        certainty   --  How sure you want to be than n really is a prime. See
                        below for details.

    Returns:

        True        --  n is probably prime. The probability that n is not
                        prime is <= 1/2**certainty.
        False       --  n is *definitely* not prime.
    """


    rand = random.SystemRandom()
    def test_once(d,r,n):
        a = rand.randint(2, n-2)
        x = pow(a, d, n)
        if x == 1 or x == n-1: return True
        for i in range(1,r):
            x = (x*x) % n
            if x == n-1: return True
        return False

    # Take care of special cases
    if n == 2 or n == 3: return True
    if n == 1 or n % 2 == 0: return False

    # Find r, d such that n = 2**r * d + 1, with d odd
    d = n - 1
    r = 0
    while d % 2 == 0:
        d = d // 2
        r += 1
    for t in range(certainty//2+1):
        if test_once(d,r,n) == False:
            return False
    return True

def find_prime(start=0, bits=0, linear=True):
    """Find a randomly chosen prime of a given number of bits, or find the
    first prime >= a given number.

    Keyword arguments:

        bits    --  Number of bits if you want a random prime.

        start   --  Starting search position if you want a non-random prime.

        linear  --  Controls search for random primes. See below.

        Note: you should specify exactly one of "bits" and "start".

    Returns:

        prime   --  A prime number

    When searching for a random prime, "linear" specifies what happens if
    the random number it picks is not prime. If linear is True, it will add
    2 and check that number. If linear is False, it will pick another random
    number to check. When linear is True, if adding 2 would take it past
    2**bits, it will wrap to 2**(bits-1)+1.
    """
    rand = random.SystemRandom()
    if bits != 0 and start == 0:
        base = 2**(bits-1)
        wrap = 2*base
        p = rand.randint(base, wrap-1) | 1
    elif bits == 0 and start != 0:
        if start <= 2:
            return 2
        p = start | 1
        linear = True
        base = start ^ 1
        wrap = 2 * p
    else:
        print("Error: find_prime must be given exactly one of bits and start")
        return 0
    while not rabin_miller(p):
        if linear:
            p += 2
            if p >= wrap:
                p = base + 1
        else:
            p = rand.randint(base, wrap-1) | 1
    return p


if __name__ == "__main__":
    cores = 1
    worker = 0
    if len(sys.argv) == 3:
        cores = int(sys.argv[1])
        worker = int(sys.argv[2])
    start = 2**1024+1
    p = find_prime(start=start)
    if worker == 0: print("p =", p)
    print(worker, rabin_miller_deterministic(p, cores=cores, worker=worker))