You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
614 lines
16 KiB
614 lines
16 KiB
#! /usr/bin/env python |
|
# |
|
# Provide some simple capabilities from number theory. |
|
# |
|
# Version of 2008.11.14. |
|
# |
|
# Written in 2005 and 2006 by Peter Pearson and placed in the public domain. |
|
# Revision history: |
|
# 2008.11.14: Use pow( base, exponent, modulus ) for modular_exp. |
|
# Make gcd and lcm accept arbitrarly many arguments. |
|
|
|
|
|
|
|
import math |
|
import types |
|
|
|
|
|
class Error( Exception ): |
|
"""Base class for exceptions in this module.""" |
|
pass |
|
|
|
class SquareRootError( Error ): |
|
pass |
|
|
|
class NegativeExponentError( Error ): |
|
pass |
|
|
|
|
|
def modular_exp( base, exponent, modulus ): |
|
"Raise base to exponent, reducing by modulus" |
|
if exponent < 0: |
|
raise NegativeExponentError( "Negative exponents (%d) not allowed" \ |
|
% exponent ) |
|
return pow( base, exponent, modulus ) |
|
# result = 1L |
|
# x = exponent |
|
# b = base + 0L |
|
# while x > 0: |
|
# if x % 2 > 0: result = (result * b) % modulus |
|
# x = x // 2 |
|
# b = ( b * b ) % modulus |
|
# return result |
|
|
|
|
|
def polynomial_reduce_mod( poly, polymod, p ): |
|
"""Reduce poly by polymod, integer arithmetic modulo p. |
|
|
|
Polynomials are represented as lists of coefficients |
|
of increasing powers of x.""" |
|
|
|
# This module has been tested only by extensive use |
|
# in calculating modular square roots. |
|
|
|
# Just to make this easy, require a monic polynomial: |
|
assert polymod[-1] == 1 |
|
|
|
assert len( polymod ) > 1 |
|
|
|
while len( poly ) >= len( polymod ): |
|
if poly[-1] != 0: |
|
for i in range( 2, len( polymod ) + 1 ): |
|
poly[-i] = ( poly[-i] - poly[-1] * polymod[-i] ) % p |
|
poly = poly[0:-1] |
|
|
|
return poly |
|
|
|
|
|
|
|
def polynomial_multiply_mod( m1, m2, polymod, p ): |
|
"""Polynomial multiplication modulo a polynomial over ints mod p. |
|
|
|
Polynomials are represented as lists of coefficients |
|
of increasing powers of x.""" |
|
|
|
# This is just a seat-of-the-pants implementation. |
|
|
|
# This module has been tested only by extensive use |
|
# in calculating modular square roots. |
|
|
|
# Initialize the product to zero: |
|
|
|
prod = ( len( m1 ) + len( m2 ) - 1 ) * [0] |
|
|
|
# Add together all the cross-terms: |
|
|
|
for i in range( len( m1 ) ): |
|
for j in range( len( m2 ) ): |
|
prod[i+j] = ( prod[i+j] + m1[i] * m2[j] ) % p |
|
|
|
return polynomial_reduce_mod( prod, polymod, p ) |
|
|
|
|
|
|
|
|
|
def polynomial_exp_mod( base, exponent, polymod, p ): |
|
"""Polynomial exponentiation modulo a polynomial over ints mod p. |
|
|
|
Polynomials are represented as lists of coefficients |
|
of increasing powers of x.""" |
|
|
|
# Based on the Handbook of Applied Cryptography, algorithm 2.227. |
|
|
|
# This module has been tested only by extensive use |
|
# in calculating modular square roots. |
|
|
|
assert exponent < p |
|
|
|
if exponent == 0: return [ 1 ] |
|
|
|
G = base |
|
k = exponent |
|
if k%2 == 1: s = G |
|
else: s = [ 1 ] |
|
|
|
while k > 1: |
|
k = k // 2 |
|
G = polynomial_multiply_mod( G, G, polymod, p ) |
|
if k%2 == 1: s = polynomial_multiply_mod( G, s, polymod, p ) |
|
|
|
return s |
|
|
|
|
|
|
|
def jacobi( a, n ): |
|
"""Jacobi symbol""" |
|
|
|
# Based on the Handbook of Applied Cryptography (HAC), algorithm 2.149. |
|
|
|
# This function has been tested by comparison with a small |
|
# table printed in HAC, and by extensive use in calculating |
|
# modular square roots. |
|
|
|
assert n >= 3 |
|
assert n%2 == 1 |
|
a = a % n |
|
if a == 0: return 0 |
|
if a == 1: return 1 |
|
a1, e = a, 0 |
|
while a1%2 == 0: |
|
a1, e = a1//2, e+1 |
|
if e%2 == 0 or n%8 == 1 or n%8 == 7: s = 1 |
|
else: s = -1 |
|
if a1 == 1: return s |
|
if n%4 == 3 and a1%4 == 3: s = -s |
|
return s * jacobi( n % a1, a1 ) |
|
|
|
|
|
|
|
|
|
def square_root_mod_prime( a, p ): |
|
"""Modular square root of a, mod p, p prime.""" |
|
|
|
# Based on the Handbook of Applied Cryptography, algorithms 3.34 to 3.39. |
|
|
|
# This module has been tested for all values in [0,p-1] for |
|
# every prime p from 3 to 1229. |
|
|
|
assert 0 <= a < p |
|
assert 1 < p |
|
|
|
if a == 0: return 0 |
|
if p == 2: return a |
|
|
|
jac = jacobi( a, p ) |
|
if jac == -1: raise SquareRootError( "%d has no square root modulo %d" \ |
|
% ( a, p ) ) |
|
|
|
if p % 4 == 3: return modular_exp( a, (p+1)//4, p ) |
|
|
|
if p % 8 == 5: |
|
d = modular_exp( a, (p-1)//4, p ) |
|
if d == 1: return modular_exp( a, (p+3)//8, p ) |
|
if d == p-1: return ( 2 * a * modular_exp( 4*a, (p-5)//8, p ) ) % p |
|
raise RuntimeError, "Shouldn't get here." |
|
|
|
for b in range( 2, p ): |
|
if jacobi( b*b-4*a, p ) == -1: |
|
f = ( a, -b, 1 ) |
|
ff = polynomial_exp_mod( ( 0, 1 ), (p+1)//2, f, p ) |
|
assert ff[1] == 0 |
|
return ff[0] |
|
raise RuntimeError, "No b found." |
|
|
|
|
|
|
|
def inverse_mod( a, m ): |
|
"""Inverse of a mod m.""" |
|
|
|
if a < 0 or m <= a: a = a % m |
|
|
|
# From Ferguson and Schneier, roughly: |
|
|
|
c, d = a, m |
|
uc, vc, ud, vd = 1, 0, 0, 1 |
|
while c != 0: |
|
q, c, d = divmod( d, c ) + ( c, ) |
|
uc, vc, ud, vd = ud - q*uc, vd - q*vc, uc, vc |
|
|
|
# At this point, d is the GCD, and ud*a+vd*m = d. |
|
# If d == 1, this means that ud is a inverse. |
|
|
|
assert d == 1 |
|
if ud > 0: return ud |
|
else: return ud + m |
|
|
|
|
|
def gcd2(a, b): |
|
"""Greatest common divisor using Euclid's algorithm.""" |
|
while a: |
|
a, b = b%a, a |
|
return b |
|
|
|
|
|
def gcd( *a ): |
|
"""Greatest common divisor. |
|
|
|
Usage: gcd( [ 2, 4, 6 ] ) |
|
or: gcd( 2, 4, 6 ) |
|
""" |
|
|
|
if len( a ) > 1: return reduce( gcd2, a ) |
|
if hasattr( a[0], "__iter__" ): return reduce( gcd2, a[0] ) |
|
return a[0] |
|
|
|
|
|
def lcm2(a,b): |
|
"""Least common multiple of two integers.""" |
|
|
|
return (a*b)//gcd(a,b) |
|
|
|
|
|
def lcm( *a ): |
|
"""Least common multiple. |
|
|
|
Usage: lcm( [ 3, 4, 5 ] ) |
|
or: lcm( 3, 4, 5 ) |
|
""" |
|
|
|
if len( a ) > 1: return reduce( lcm2, a ) |
|
if hasattr( a[0], "__iter__" ): return reduce( lcm2, a[0] ) |
|
return a[0] |
|
|
|
|
|
|
|
def factorization( n ): |
|
"""Decompose n into a list of (prime,exponent) pairs.""" |
|
|
|
assert isinstance( n, types.IntType ) or isinstance( n, types.LongType ) |
|
|
|
if n < 2: return [] |
|
|
|
result = [] |
|
d = 2 |
|
|
|
# Test the small primes: |
|
|
|
for d in smallprimes: |
|
if d > n: break |
|
q, r = divmod( n, d ) |
|
if r == 0: |
|
count = 1 |
|
while d <= n: |
|
n = q |
|
q, r = divmod( n, d ) |
|
if r != 0: break |
|
count = count + 1 |
|
result.append( ( d, count ) ) |
|
|
|
# If n is still greater than the last of our small primes, |
|
# it may require further work: |
|
|
|
if n > smallprimes[-1]: |
|
if is_prime( n ): # If what's left is prime, it's easy: |
|
result.append( ( n, 1 ) ) |
|
else: # Ugh. Search stupidly for a divisor: |
|
d = smallprimes[-1] |
|
while 1: |
|
d = d + 2 # Try the next divisor. |
|
q, r = divmod( n, d ) |
|
if q < d: break # n < d*d means we're done, n = 1 or prime. |
|
if r == 0: # d divides n. How many times? |
|
count = 1 |
|
n = q |
|
while d <= n: # As long as d might still divide n, |
|
q, r = divmod( n, d ) # see if it does. |
|
if r != 0: break |
|
n = q # It does. Reduce n, increase count. |
|
count = count + 1 |
|
result.append( ( d, count ) ) |
|
if n > 1: result.append( ( n, 1 ) ) |
|
|
|
return result |
|
|
|
|
|
|
|
def phi( n ): |
|
"""Return the Euler totient function of n.""" |
|
|
|
assert isinstance( n, types.IntType ) or isinstance( n, types.LongType ) |
|
|
|
if n < 3: return 1 |
|
|
|
result = 1 |
|
ff = factorization( n ) |
|
for f in ff: |
|
e = f[1] |
|
if e > 1: |
|
result = result * f[0] ** (e-1) * ( f[0] - 1 ) |
|
else: |
|
result = result * ( f[0] - 1 ) |
|
return result |
|
|
|
|
|
def carmichael( n ): |
|
"""Return Carmichael function of n. |
|
|
|
Carmichael(n) is the smallest integer x such that |
|
m**x = 1 mod n for all m relatively prime to n. |
|
""" |
|
|
|
return carmichael_of_factorized( factorization( n ) ) |
|
|
|
|
|
def carmichael_of_factorized( f_list ): |
|
"""Return the Carmichael function of a number that is |
|
represented as a list of (prime,exponent) pairs. |
|
""" |
|
|
|
if len( f_list ) < 1: return 1 |
|
|
|
result = carmichael_of_ppower( f_list[0] ) |
|
for i in range( 1, len( f_list ) ): |
|
result = lcm( result, carmichael_of_ppower( f_list[i] ) ) |
|
|
|
return result |
|
|
|
def carmichael_of_ppower( pp ): |
|
"""Carmichael function of the given power of the given prime. |
|
""" |
|
|
|
p, a = pp |
|
if p == 2 and a > 2: return 2**(a-2) |
|
else: return (p-1) * p**(a-1) |
|
|
|
|
|
|
|
def order_mod( x, m ): |
|
"""Return the order of x in the multiplicative group mod m. |
|
""" |
|
|
|
# Warning: this implementation is not very clever, and will |
|
# take a long time if m is very large. |
|
|
|
if m <= 1: return 0 |
|
|
|
assert gcd( x, m ) == 1 |
|
|
|
z = x |
|
result = 1 |
|
while z != 1: |
|
z = ( z * x ) % m |
|
result = result + 1 |
|
return result |
|
|
|
|
|
def largest_factor_relatively_prime( a, b ): |
|
"""Return the largest factor of a relatively prime to b. |
|
""" |
|
|
|
while 1: |
|
d = gcd( a, b ) |
|
if d <= 1: break |
|
b = d |
|
while 1: |
|
q, r = divmod( a, d ) |
|
if r > 0: |
|
break |
|
a = q |
|
return a |
|
|
|
|
|
def kinda_order_mod( x, m ): |
|
"""Return the order of x in the multiplicative group mod m', |
|
where m' is the largest factor of m relatively prime to x. |
|
""" |
|
|
|
return order_mod( x, largest_factor_relatively_prime( m, x ) ) |
|
|
|
|
|
def is_prime( n ): |
|
"""Return True if x is prime, False otherwise. |
|
|
|
We use the Miller-Rabin test, as given in Menezes et al. p. 138. |
|
This test is not exact: there are composite values n for which |
|
it returns True. |
|
|
|
In testing the odd numbers from 10000001 to 19999999, |
|
about 66 composites got past the first test, |
|
5 got past the second test, and none got past the third. |
|
Since factors of 2, 3, 5, 7, and 11 were detected during |
|
preliminary screening, the number of numbers tested by |
|
Miller-Rabin was (19999999 - 10000001)*(2/3)*(4/5)*(6/7) |
|
= 4.57 million. |
|
""" |
|
|
|
# (This is used to study the risk of false positives:) |
|
global miller_rabin_test_count |
|
|
|
miller_rabin_test_count = 0 |
|
|
|
if n <= smallprimes[-1]: |
|
if n in smallprimes: return True |
|
else: return False |
|
|
|
if gcd( n, 2*3*5*7*11 ) != 1: return False |
|
|
|
# Choose a number of iterations sufficient to reduce the |
|
# probability of accepting a composite below 2**-80 |
|
# (from Menezes et al. Table 4.4): |
|
|
|
t = 40 |
|
n_bits = 1 + int( math.log( n, 2 ) ) |
|
for k, tt in ( ( 100, 27 ), |
|
( 150, 18 ), |
|
( 200, 15 ), |
|
( 250, 12 ), |
|
( 300, 9 ), |
|
( 350, 8 ), |
|
( 400, 7 ), |
|
( 450, 6 ), |
|
( 550, 5 ), |
|
( 650, 4 ), |
|
( 850, 3 ), |
|
( 1300, 2 ), |
|
): |
|
if n_bits < k: break |
|
t = tt |
|
|
|
# Run the test t times: |
|
|
|
s = 0 |
|
r = n - 1 |
|
while ( r % 2 ) == 0: |
|
s = s + 1 |
|
r = r // 2 |
|
for i in xrange( t ): |
|
a = smallprimes[ i ] |
|
y = modular_exp( a, r, n ) |
|
if y != 1 and y != n-1: |
|
j = 1 |
|
while j <= s - 1 and y != n - 1: |
|
y = modular_exp( y, 2, n ) |
|
if y == 1: |
|
miller_rabin_test_count = i + 1 |
|
return False |
|
j = j + 1 |
|
if y != n-1: |
|
miller_rabin_test_count = i + 1 |
|
return False |
|
return True |
|
|
|
|
|
def next_prime( starting_value ): |
|
"Return the smallest prime larger than the starting value." |
|
|
|
if starting_value < 2: return 2 |
|
result = ( starting_value + 1 ) | 1 |
|
while not is_prime( result ): result = result + 2 |
|
return result |
|
|
|
|
|
smallprimes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, |
|
43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, |
|
101, 103, 107, 109, 113, 127, 131, 137, 139, 149, |
|
151, 157, 163, 167, 173, 179, 181, 191, 193, 197, |
|
199, 211, 223, 227, 229, 233, 239, 241, 251, 257, |
|
263, 269, 271, 277, 281, 283, 293, 307, 311, 313, |
|
317, 331, 337, 347, 349, 353, 359, 367, 373, 379, |
|
383, 389, 397, 401, 409, 419, 421, 431, 433, 439, |
|
443, 449, 457, 461, 463, 467, 479, 487, 491, 499, |
|
503, 509, 521, 523, 541, 547, 557, 563, 569, 571, |
|
577, 587, 593, 599, 601, 607, 613, 617, 619, 631, |
|
641, 643, 647, 653, 659, 661, 673, 677, 683, 691, |
|
701, 709, 719, 727, 733, 739, 743, 751, 757, 761, |
|
769, 773, 787, 797, 809, 811, 821, 823, 827, 829, |
|
839, 853, 857, 859, 863, 877, 881, 883, 887, 907, |
|
911, 919, 929, 937, 941, 947, 953, 967, 971, 977, |
|
983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, |
|
1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, |
|
1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, |
|
1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229] |
|
|
|
miller_rabin_test_count = 0 |
|
|
|
def __main__(): |
|
|
|
# Making sure locally defined exceptions work: |
|
# p = modular_exp( 2, -2, 3 ) |
|
# p = square_root_mod_prime( 2, 3 ) |
|
|
|
|
|
print "Testing gcd..." |
|
assert gcd( 3*5*7, 3*5*11, 3*5*13 ) == 3*5 |
|
assert gcd( [ 3*5*7, 3*5*11, 3*5*13 ] ) == 3*5 |
|
assert gcd( 3 ) == 3 |
|
|
|
print "Testing lcm..." |
|
assert lcm( 3, 5*3, 7*3 ) == 3*5*7 |
|
assert lcm( [ 3, 5*3, 7*3 ] ) == 3*5*7 |
|
assert lcm( 3 ) == 3 |
|
|
|
print "Testing next_prime..." |
|
bigprimes = ( 999671, |
|
999683, |
|
999721, |
|
999727, |
|
999749, |
|
999763, |
|
999769, |
|
999773, |
|
999809, |
|
999853, |
|
999863, |
|
999883, |
|
999907, |
|
999917, |
|
999931, |
|
999953, |
|
999959, |
|
999961, |
|
999979, |
|
999983 ) |
|
|
|
for i in xrange( len( bigprimes ) - 1 ): |
|
assert next_prime( bigprimes[i] ) == bigprimes[ i+1 ] |
|
|
|
error_tally = 0 |
|
|
|
# Test the square_root_mod_prime function: |
|
|
|
for p in smallprimes: |
|
print "Testing square_root_mod_prime for modulus p = %d." % p |
|
squares = [] |
|
|
|
for root in range( 0, 1+p//2 ): |
|
sq = ( root * root ) % p |
|
squares.append( sq ) |
|
calculated = square_root_mod_prime( sq, p ) |
|
if ( calculated * calculated ) % p != sq: |
|
error_tally = error_tally + 1 |
|
print "Failed to find %d as sqrt( %d ) mod %d. Said %d." % \ |
|
( root, sq, p, calculated ) |
|
|
|
for nonsquare in range( 0, p ): |
|
if nonsquare not in squares: |
|
try: |
|
calculated = square_root_mod_prime( nonsquare, p ) |
|
except SquareRootError: |
|
pass |
|
else: |
|
error_tally = error_tally + 1 |
|
print "Failed to report no root for sqrt( %d ) mod %d." % \ |
|
( nonsquare, p ) |
|
|
|
# Test the jacobi function: |
|
for m in range( 3, 400, 2 ): |
|
print "Testing jacobi for modulus m = %d." % m |
|
if is_prime( m ): |
|
squares = [] |
|
for root in range( 1, m ): |
|
if jacobi( root * root, m ) != 1: |
|
error_tally = error_tally + 1 |
|
print "jacobi( %d * %d, %d ) != 1" % ( root, root, m ) |
|
squares.append( root * root % m ) |
|
for i in range( 1, m ): |
|
if not i in squares: |
|
if jacobi( i, m ) != -1: |
|
error_tally = error_tally + 1 |
|
print "jacobi( %d, %d ) != -1" % ( i, m ) |
|
else: # m is not prime. |
|
f = factorization( m ) |
|
for a in range( 1, m ): |
|
c = 1 |
|
for i in f: |
|
c = c * jacobi( a, i[0] ) ** i[1] |
|
if c != jacobi( a, m ): |
|
error_tally = error_tally + 1 |
|
print "%d != jacobi( %d, %d )" % ( c, a, m ) |
|
|
|
|
|
# Test the inverse_mod function: |
|
print "Testing inverse_mod . . ." |
|
import random |
|
n_tests = 0 |
|
for i in range( 100 ): |
|
m = random.randint( 20, 10000 ) |
|
for j in range( 100 ): |
|
a = random.randint( 1, m-1 ) |
|
if gcd( a, m ) == 1: |
|
n_tests = n_tests + 1 |
|
inv = inverse_mod( a, m ) |
|
if inv <= 0 or inv >= m or ( a * inv ) % m != 1: |
|
error_tally = error_tally + 1 |
|
print "%d = inverse_mod( %d, %d ) is wrong." % ( inv, a, m ) |
|
assert n_tests > 1000 |
|
print n_tests, " tests of inverse_mod completed." |
|
|
|
class FailedTest(Exception): pass |
|
print error_tally, "errors detected." |
|
if error_tally != 0: |
|
raise FailedTest("%d errors detected" % error_tally) |
|
|
|
if __name__ == '__main__': |
|
__main__()
|
|
|