Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- """
- 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
- """
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement