1 | # -*- coding: utf-8 -*- |
||
2 | # |
||
3 | # Copyright 2011 Sybren A. Stüvel <[email protected]> |
||
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 | View Code Duplication | def jacobi(a, b): |
|
0 ignored issues
–
show
Duplication
introduced
by
![]() |
|||
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 | View Code Duplication | def jacobi_witness(x, n): |
|
0 ignored issues
–
show
|
|||
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 | View Code Duplication | if __name__ == '__main__': |
|
0 ignored issues
–
show
|
|||
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' |
||
149 |