#### 167 lines 3.8 KiB Python Raw Blame History

 ```# -*- coding: utf-8 -*- ``` ```# ``` ```# Copyright 2011 Sybren A. Stüvel ``` ```# ``` ```# 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') ```