 # Frobenius primality test

Nov 13th, 2012
348
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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.    """