rsa_source/rsa/prime.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 #  Copyright 2011 Sybren A. Stüvel <sybren@stuvel.eu>
   4 #
   5 #  Licensed under the Apache License, Version 2.0 (the "License");
   6 #  you may not use this file except in compliance with the License.
   7 #  You may obtain a copy of the License at
   8 #
   9 #      http://www.apache.org/licenses/LICENSE-2.0
  10 #
  11 #  Unless required by applicable law or agreed to in writing, software
  12 #  distributed under the License is distributed on an "AS IS" BASIS,
  13 #  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14 #  See the License for the specific language governing permissions and
  15 #  limitations under the License.
  16
  17 '''Numerical functions related to primes.'''
  18
  19 __all__ = [ 'getprime', 'are_relatively_prime']
  20
  21 import rsa.randnum
  22
  23 def gcd(p, q):
  24     """Returns the greatest common divisor of p and q
  25
  26     >>> gcd(48, 180)
  27     12
  28     """
  29
  30     while q != 0:
  31         if p < q: (p,q) = (q,p)
  32         (p,q) = (q, p % q)
  33     return p
  34
  35
  36 def jacobi(a, b):
  37     """Calculates the value of the Jacobi symbol (a/b) where both a and b are
  38     positive integers, and b is odd
  39     """
  40
  41     if a == 0: return 0
  42     result = 1
  43     while a > 1:
  44         if a & 1:
  45             if ((a-1)*(b-1) >> 2) & 1:
  46                 result = -result
  47             a, b = b % a, a
  48         else:
  49             if (((b * b) - 1) >> 3) & 1:
  50                 result = -result
  51             a >>= 1
  52     if a == 0: return 0
  53     return result
  54
  55 def jacobi_witness(x, n):
  56     """Returns False if n is an Euler pseudo-prime with base x, and
  57     True otherwise.
  58     """
  59
  60     j = jacobi(x, n) % n
  61     f = pow(x, (n - 1) // 2, n)
  62
  63     if j == f: return False
  64     return True
  65
  66 def randomized_primality_testing(n, k):
  67     """Calculates whether n is composite (which is always correct) or
  68     prime (which is incorrect with error probability 2**-k)
  69
  70     Returns False if the number is composite, and True if it's
  71     probably prime.
  72     """
  73
  74     # 50% of Jacobi-witnesses can report compositness of non-prime numbers
  75
  76     for _ in range(k):
  77         x = rsa.randnum.randint(n-1)
  78         if jacobi_witness(x, n): return False
  79
  80     return True
  81
  82 def is_prime(number):
  83     """Returns True if the number is prime, and False otherwise.
  84
  85     >>> is_prime(42)
  86     False
  87     >>> is_prime(41)
  88     True
  89     """
  90
  91     return randomized_primality_testing(number, 6)
  92
  93 def getprime(nbits):
  94     """Returns a prime number that can be stored in 'nbits' bits.
  95
  96     >>> p = getprime(128)
  97     >>> is_prime(p-1)
  98     False
  99     >>> is_prime(p)
 100     True
 101     >>> is_prime(p+1)
 102     False
 103
 104     >>> from rsa import common
 105     >>> common.bit_size(p) <= 128
 106     True
 107
 108     """
 109
 110     while True:
 111         integer = rsa.randnum.read_random_int(nbits)
 112
 113         # Make sure it's odd
 114         integer |= 1
 115
 116         # Test for primeness
 117         if is_prime(integer):
 118             return integer
 119
 120         # Retry if not prime
 121
 122
 123 def are_relatively_prime(a, b):
 124     """Returns True if a and b are relatively prime, and False if they
 125     are not.
 126
 127     >>> are_relatively_prime(2, 3)
 128     1
 129     >>> are_relatively_prime(2, 4)
 130     0
 131     """
 132
 133     d = gcd(a, b)
 134     return (d == 1)
 135
 136 if __name__ == '__main__':
 137     print 'Running doctests 1000x or until failure'
 138     import doctest
 139
 140     for count in range(1000):
 141         (failures, tests) = doctest.testmod()
 142         if failures:
 143             break
 144
 145         if count and count % 100 == 0:
 146             print '%i times' % count
 147
 148     print 'Doctests done'