PyBitmessage/rsa/prime.py

167 lines
3.8 KiB
Python
Raw Normal View History

2012-11-19 20:45:05 +01:00
# -*- 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')