167 lines
3.8 KiB
Python
167 lines
3.8 KiB
Python
# -*- coding: utf-8 -*-
|
|
#
|
|
# Copyright 2011 Sybren A. Stüvel <sybren@stuvel.eu>
|
|
#
|
|
# Licensed under the Apache License, Version 2.0 (the "License");
|
|
# you may not use this file except in compliance with the License.
|
|
# You may obtain a copy of the License at
|
|
#
|
|
# http://www.apache.org/licenses/LICENSE-2.0
|
|
#
|
|
# Unless required by applicable law or agreed to in writing, software
|
|
# distributed under the License is distributed on an "AS IS" BASIS,
|
|
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
# See the License for the specific language governing permissions and
|
|
# limitations under the License.
|
|
|
|
'''Numerical functions related to primes.
|
|
|
|
Implementation based on the book Algorithm Design by Michael T. Goodrich and
|
|
Roberto Tamassia, 2002.
|
|
'''
|
|
|
|
__all__ = [ 'getprime', 'are_relatively_prime']
|
|
|
|
import rsa.randnum
|
|
|
|
def gcd(p, q):
|
|
'''Returns the greatest common divisor of p and q
|
|
|
|
>>> gcd(48, 180)
|
|
12
|
|
'''
|
|
|
|
while q != 0:
|
|
if p < q: (p,q) = (q,p)
|
|
(p,q) = (q, p % q)
|
|
return p
|
|
|
|
|
|
def jacobi(a, b):
|
|
'''Calculates the value of the Jacobi symbol (a/b) where both a and b are
|
|
positive integers, and b is odd
|
|
|
|
:returns: -1, 0 or 1
|
|
'''
|
|
|
|
assert a > 0
|
|
assert b > 0
|
|
|
|
if a == 0: return 0
|
|
result = 1
|
|
while a > 1:
|
|
if a & 1:
|
|
if ((a-1)*(b-1) >> 2) & 1:
|
|
result = -result
|
|
a, b = b % a, a
|
|
else:
|
|
if (((b * b) - 1) >> 3) & 1:
|
|
result = -result
|
|
a >>= 1
|
|
if a == 0: return 0
|
|
return result
|
|
|
|
def jacobi_witness(x, n):
|
|
'''Returns False if n is an Euler pseudo-prime with base x, and
|
|
True otherwise.
|
|
'''
|
|
|
|
j = jacobi(x, n) % n
|
|
|
|
f = pow(x, n >> 1, n)
|
|
|
|
if j == f: return False
|
|
return True
|
|
|
|
def randomized_primality_testing(n, k):
|
|
'''Calculates whether n is composite (which is always correct) or
|
|
prime (which is incorrect with error probability 2**-k)
|
|
|
|
Returns False if the number is composite, and True if it's
|
|
probably prime.
|
|
'''
|
|
|
|
# 50% of Jacobi-witnesses can report compositness of non-prime numbers
|
|
|
|
# The implemented algorithm using the Jacobi witness function has error
|
|
# probability q <= 0.5, according to Goodrich et. al
|
|
#
|
|
# q = 0.5
|
|
# t = int(math.ceil(k / log(1 / q, 2)))
|
|
# So t = k / log(2, 2) = k / 1 = k
|
|
# this means we can use range(k) rather than range(t)
|
|
|
|
for _ in range(k):
|
|
x = rsa.randnum.randint(n-1)
|
|
if jacobi_witness(x, n): return False
|
|
|
|
return True
|
|
|
|
def is_prime(number):
|
|
'''Returns True if the number is prime, and False otherwise.
|
|
|
|
>>> is_prime(42)
|
|
False
|
|
>>> is_prime(41)
|
|
True
|
|
'''
|
|
|
|
return randomized_primality_testing(number, 6)
|
|
|
|
def getprime(nbits):
|
|
'''Returns a prime number that can be stored in 'nbits' bits.
|
|
|
|
>>> p = getprime(128)
|
|
>>> is_prime(p-1)
|
|
False
|
|
>>> is_prime(p)
|
|
True
|
|
>>> is_prime(p+1)
|
|
False
|
|
|
|
>>> from rsa import common
|
|
>>> common.bit_size(p) == 128
|
|
True
|
|
|
|
'''
|
|
|
|
while True:
|
|
integer = rsa.randnum.read_random_int(nbits)
|
|
|
|
# Make sure it's odd
|
|
integer |= 1
|
|
|
|
# Test for primeness
|
|
if is_prime(integer):
|
|
return integer
|
|
|
|
# Retry if not prime
|
|
|
|
|
|
def are_relatively_prime(a, b):
|
|
'''Returns True if a and b are relatively prime, and False if they
|
|
are not.
|
|
|
|
>>> are_relatively_prime(2, 3)
|
|
1
|
|
>>> are_relatively_prime(2, 4)
|
|
0
|
|
'''
|
|
|
|
d = gcd(a, b)
|
|
return (d == 1)
|
|
|
|
if __name__ == '__main__':
|
|
print('Running doctests 1000x or until failure')
|
|
import doctest
|
|
|
|
for count in range(1000):
|
|
(failures, tests) = doctest.testmod()
|
|
if failures:
|
|
break
|
|
|
|
if count and count % 100 == 0:
|
|
print('%i times' % count)
|
|
|
|
print('Doctests done')
|