1. """
  2.    Frobenius Pseudoprimality Test: Given an integer n greater than 2, determine
  3.    with high probability if it is COMPOSITE or PROBABLY PRIME:
  4.  
  5.    1. [Choose parameters] Set a and b to distinct random numbers uniformly
  6.    distributed on the range 1 to n-1 inclusive. Set d = a^2 - 4b. If d is a
  7.    perfect square or if gcd(2abd, n) 1 go to Step 1.
  8.  
  9.    2. [Initialization] Set m = (n - (d/n)) / 2, where (d/n) is the jacobi
  10.    symbol. Set w = a^2 * b^(-1) - 2 (mod n), where b^(-1) is the modular
  11.    inverse of b with respect to n. Set u = 2. Set v = 2. Set ms to a list of
  12.    the bits in the binary representation of m, most significant bit first.
  13.  
  14.    3. [Compute Lucas chain] If ms is empty, go to Step 4. If the first element
  15.    of ms is 1, set u = u * v - w (mod n), set v = v * v - 2 (mod n), remove
  16.    the first element of ms, and go to Step 3. If the first element of ms is 0,
  17.    set v = u * v - w (mod n), set u = u * u - 2 (mod n), remove the first
  18.    element of ms, and go to Step 3.
  19.  
  20.    4. [Lucas test] If w * u 2 * v (mod n), report COMPOSITE and halt.
  21.  
  22.    5. [Frobenius test] Set x = b ^ ((n-1)/2) (mod n). If x * u == 2 (mod n),
  23.    report PROBABLY PRIME, otherwise report COMPOSITE. Halt.
  24. """
  25. from random import randint
  26. from fractions import gcd
  27. import math
  28. from ma.primec import inverse, bits_of
  29. import sys
  30.  
  31. def jacobi(a,b):
  32.     if b <= 0 or b % 2 ==0:
  33.         return 0
  34.     j = 1
  35.     if a < 0:
  36.         a = -a
  37.         if (b % 4) == 3:
  38.             j = -j
  39.     while a != 0:
  40.         while (a % 2) == 0:
  41.             # Process factors of 2: Jacobi(2,b)=-1 if b=3,5 (mod 8)
  42.             a //= 2
  43.             if (b % 8) in [3, 5]:
  44.                 j = -j
  45.         # Quadratic reciprocity: Jacobi(a,b)=-Jacobi(b,a) if a=3,b=3 (mod 4)
  46.         a, b = b, a
  47.         if (a % 4)==3 and (b % 4)==3:
  48.             j = -j
  49.         a = a % b
  50.     return j if b == 1 else 0
  51.  
  52. def is_square(n):
  53.     if n < 0:
  54.         return False
  55.     return int(math.sqrt(n+0.5)) ** 2 == n
  56.  
  57. def frobenius(n):
  58.     """ returns True if n is probably prime
  59.    """
  60.     if n % 2 == 0 or n < 2:
  61.         return False
  62.     #step 1
  63.     while 1:
  64.         a = randint(1, n-1)
  65.         b = randint(1, n-1)
  66.         d = a ** 2 - 4 * b
  67.         if  not (is_square(d) or gcd(2 * a  * b * d, n) != 1):
  68.             break
  69.     #step 2
  70.     m = (n - jacobi(d, n)) // 2
  71.     v = w = (a ** 2  * inverse(b, n) - 2) % n
  72.     u = 2
  73.     # step 3
  74.     for bit in bits_of(m):
  75.         if bit:
  76.             u = (u * v - w) % n
  77.             v = (v * v - 2) % n
  78.         else:
  79.             v = (u * v - w) % n
  80.             u = (u * u - 2) % n
  81.     # step 4
  82.     if (w * u - 2 * v) % n != 0:
  83.         return False
  84.     return (pow(b, (n - 1) // 2, n) * u) % n ==  2
  85.        
  86. def test_jacobi():
  87.     cases = [(1101, 9907, -1), (1001, 9907, -1), (1236, 20003, 1),
  88.         (1235, 20003, -1), (37603, 48611, 1), (40, 17, -1) ,(40, 11, -1),
  89.         (0, 3, 0)]
  90.     error = [ (0, 1, 0)]
  91.     for m, n, e in cases:
  92.         assert jacobi(m,n) == e
  93.  
  94. def test_inv():
  95.     cases = [(29, 35, 78), (7,8,11)]
  96.     for a, i, m in cases:
  97.         assert inverse(a, m) == i
  98.         assert inverse(i, m) == a
  99.  
  100. def test_sq():
  101.     squres = [i**2 for i in range(10)]
  102.     def sq(n):
  103.         return n in squres
  104.     for i in range(100):
  105.         assert sq(i) == is_square(i)
  106.        
  107.            
  108. if __name__ == '__main__':
  109.     for i in range(3, 100, 2):
  110.         if frobenius(i):
  111.             print i,
  112.     print
  113.     print frobenius(3825123056546413051)
  114.     print frobenius(2**89 - 1)
  115.     """ 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97
  116.        False
  117.        True
  118.    """