PyBitmessage/rsa/prime.py

# -*- 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')