"""
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
"""