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