/ adafruit_rsa / prime.py
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  #      https://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  Implementation based on the book Algorithm Design by Michael T. Goodrich and
 20  Roberto Tamassia, 2002.
 21  """
 22  # pylint: disable=invalid-name
 23  import adafruit_rsa.common
 24  import adafruit_rsa.randnum
 25  
 26  __version__ = "0.0.0-auto.0"
 27  __repo__ = "https://github.com/adafruit/Adafruit_CircuitPython_RSA.git"
 28  
 29  __all__ = ["getprime", "are_relatively_prime"]
 30  
 31  
 32  def gcd(p, q):
 33      """Returns the greatest common divisor of p and q
 34  
 35      >>> gcd(48, 180)
 36      12
 37      """
 38  
 39      while q != 0:
 40          (p, q) = (q, p % q)
 41      return p
 42  
 43  
 44  def get_primality_testing_rounds(number):
 45      """Returns minimum number of rounds for Miller-Rabing primality testing,
 46      based on number bitsize.
 47  
 48      According to NIST FIPS 186-4, Appendix C, Table C.3, minimum number of
 49      rounds of M-R testing, using an error probability of 2 ** (-100), for
 50      different p, q bitsizes are:
 51        * p, q bitsize: 512; rounds: 7
 52        * p, q bitsize: 1024; rounds: 4
 53        * p, q bitsize: 1536; rounds: 3
 54      See: http://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf
 55      """
 56  
 57      # Calculate number bitsize.
 58      bitsize = adafruit_rsa.common.bit_size(number)
 59      # Set number of rounds.
 60      if bitsize >= 1536:
 61          return 3
 62      if bitsize >= 1024:
 63          return 4
 64      if bitsize >= 512:
 65          return 7
 66      # For smaller bitsizes, set arbitrary number of rounds.
 67      return 10
 68  
 69  
 70  def miller_rabin_primality_testing(n, k):
 71      """Calculates whether n is composite (which is always correct) or prime
 72      (which theoretically is incorrect with error probability 4**-k), by
 73      applying Miller-Rabin primality testing.
 74  
 75      For reference and implementation example, see:
 76      https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
 77  
 78      :param n: Integer to be tested for primality.
 79      :type n: int
 80      :param k: Number of rounds (witnesses) of Miller-Rabin testing.
 81      :type k: int
 82      :return: False if the number is composite, True if it's probably prime.
 83      :rtype: bool
 84      """
 85  
 86      # prevent potential infinite loop when d = 0
 87      if n < 2:
 88          return False
 89  
 90      # Decompose (n - 1) to write it as (2 ** r) * d
 91      # While d is even, divide it by 2 and increase the exponent.
 92      d = n - 1
 93      r = 0
 94  
 95      while not d & 1:
 96          r += 1
 97          d >>= 1
 98  
 99      # Test k witnesses.
100      for _ in range(k):
101          # Generate random integer a, where 2 <= a <= (n - 2)
102          a = adafruit_rsa.randnum.randint(n - 3) + 1
103  
104          x = adafruit_rsa.core.fast_pow(a, d, n)
105  
106          if x in (1, n - 1):
107              continue
108  
109          for _ in range(r - 1):
110              x = adafruit_rsa.core.fast_pow(x, 2, n)
111              if x == 1:
112                  # n is composite.
113                  return False
114              if x == n - 1:
115                  # Exit inner loop and continue with next witness.
116                  break
117          else:
118              # If loop doesn't break, n is composite.
119              return False
120      return True
121  
122  
123  def pow_mod(x, y, z):
124      "Calculate (x ** y) % z efficiently."
125      number = 1
126      while y:
127          if y & 1:
128              number = number * x % z
129          y >>= 1
130          x = x * x % z
131      return number
132  
133  
134  def is_prime(number):
135      """Returns True if the number is prime, and False otherwise.
136  
137      >>> is_prime(2)
138      True
139      >>> is_prime(42)
140      False
141      >>> is_prime(41)
142      True
143      """
144  
145      # Check for small numbers.
146      if number < 10:
147          return number in {2, 3, 5, 7}
148  
149      # Check for even numbers.
150      if not number & 1:
151          return False
152  
153      # Calculate minimum number of rounds.
154      k = get_primality_testing_rounds(number)
155  
156      # Run primality testing with (minimum + 1) rounds.
157      return miller_rabin_primality_testing(number, k + 1)
158  
159  
160  def getprime(nbits):
161      """Returns a prime number that can be stored in 'nbits' bits.
162  
163      >>> p = getprime(128)
164      >>> is_prime(p-1)
165      False
166      >>> is_prime(p)
167      True
168      >>> is_prime(p+1)
169      False
170  
171      >>> from adafruit_rsa.rsa import common
172      >>> common.bit_size(p) == 128
173      True
174      """
175  
176      assert nbits > 3  # the loop wil hang on too small numbers
177  
178      while True:
179          integer = adafruit_rsa.randnum.read_random_odd_int(nbits)
180  
181          # Test for primeness
182          if is_prime(integer):
183              return integer
184  
185              # Retry if not prime
186  
187  
188  def are_relatively_prime(a, b):
189      """Returns True if a and b are relatively prime, and False if they
190      are not.
191  
192      >>> are_relatively_prime(2, 3)
193      True
194      >>> are_relatively_prime(2, 4)
195      False
196      """
197  
198      d = gcd(a, b)
199      return d == 1