"""
    Frobenius Pseudoprimality Test: Given an integer n greater than 2, determine
    with high probability if it is COMPOSITE or PROBABLY PRIME:

    1. [Choose parameters] Set a and b to distinct random numbers uniformly
    distributed on the range 1 to n-1 inclusive. Set d = a^2 - 4b. If d is a
    perfect square or if gcd(2abd, n) 1 go to Step 1.

    2. [Initialization] Set m = (n - (d/n)) / 2, where (d/n) is the jacobi
    symbol. Set w = a^2 * b^(-1) - 2 (mod n), where b^(-1) is the modular
    inverse of b with respect to n. Set u = 2. Set v = 2. Set ms to a list of
    the bits in the binary representation of m, most significant bit first.

    3. [Compute Lucas chain] If ms is empty, go to Step 4. If the first element
    of ms is 1, set u = u * v - w (mod n), set v = v * v - 2 (mod n), remove
    the first element of ms, and go to Step 3. If the first element of ms is 0,
    set v = u * v - w (mod n), set u = u * u - 2 (mod n), remove the first
    element of ms, and go to Step 3.

    4. [Lucas test] If w * u 2 * v (mod n), report COMPOSITE and halt.

    5. [Frobenius test] Set x = b ^ ((n-1)/2) (mod n). If x * u == 2 (mod n),
    report PROBABLY PRIME, otherwise report COMPOSITE. Halt.
"""
from random import randint
from fractions import gcd
import math
from ma.primec import inverse, bits_of
import sys

def jacobi(a,b):
    if b <= 0 or b % 2 ==0:
        return 0
    j = 1
    if a < 0:
        a = -a
        if (b % 4) == 3:
            j = -j
    while a != 0:
        while (a % 2) == 0:
            # Process factors of 2: Jacobi(2,b)=-1 if b=3,5 (mod 8) 
            a //= 2
            if (b % 8) in [3, 5]:
                j = -j
        # Quadratic reciprocity: Jacobi(a,b)=-Jacobi(b,a) if a=3,b=3 (mod 4) 
        a, b = b, a
        if (a % 4)==3 and (b % 4)==3:
            j = -j
        a = a % b
    return j if b == 1 else 0

def is_square(n):
    if n < 0:
        return False
    return int(math.sqrt(n+0.5)) ** 2 == n

def frobenius(n):
    """ returns True if n is probably prime
    """
    if n % 2 == 0 or n < 2:
        return False
    #step 1
    while 1:
        a = randint(1, n-1)
        b = randint(1, n-1)
        d = a ** 2 - 4 * b
        if  not (is_square(d) or gcd(2 * a  * b * d, n) != 1):
            break
    #step 2
    m = (n - jacobi(d, n)) // 2
    v = w = (a ** 2  * inverse(b, n) - 2) % n
    u = 2
    # step 3
    for bit in bits_of(m):
        if bit:
            u = (u * v - w) % n 
            v = (v * v - 2) % n
        else:
            v = (u * v - w) % n
            u = (u * u - 2) % n
    # step 4
    if (w * u - 2 * v) % n != 0:
        return False
    return (pow(b, (n - 1) // 2, n) * u) % n ==  2
        
def test_jacobi():
    cases = [(1101, 9907, -1), (1001, 9907, -1), (1236, 20003, 1),
        (1235, 20003, -1), (37603, 48611, 1), (40, 17, -1) ,(40, 11, -1),
        (0, 3, 0)]
    error = [ (0, 1, 0)]
    for m, n, e in cases:
        assert jacobi(m,n) == e
 
def test_inv():
    cases = [(29, 35, 78), (7,8,11)]
    for a, i, m in cases:
        assert inverse(a, m) == i
        assert inverse(i, m) == a

def test_sq():
    squres = [i**2 for i in range(10)]
    def sq(n):
        return n in squres
    for i in range(100):
        assert sq(i) == is_square(i)
        
            
if __name__ == '__main__':
    for i in range(3, 100, 2):
        if frobenius(i):
            print i,
    print
    print frobenius(3825123056546413051)
    print frobenius(2**89 - 1)
    """ 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
        False
        True
    """