#! /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.
from __future__ import division
from six import integer_types, PY3
from six.moves import reduce
try:
xrange
except NameError:
xrange = range
try:
from gmpy2 import powmod
GMPY2 = True
GMPY = False
except ImportError:
GMPY2 = False
try:
from gmpy import mpz
GMPY = True
except ImportError:
GMPY = False
import math
import warnings
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): # pragma: no cover
"""Raise base to exponent, reducing by modulus"""
# deprecated in 0.14
warnings.warn("Function is unused in library code. If you use this code, "
"change to pow() builtin.", DeprecationWarning)
if exponent < 0:
raise NegativeExponentError("Negative exponents (%d) not allowed"
% exponent)
return pow(base, exponent, modulus)
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 xrange(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 xrange(len(m1)):
for j in xrange(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 pow(a, (p + 1) // 4, p)
if p % 8 == 5:
d = pow(a, (p - 1) // 4, p)
if d == 1:
return pow(a, (p + 3) // 8, p)
if d == p - 1:
return (2 * a * pow(4 * a, (p - 5) // 8, p)) % p
raise RuntimeError("Shouldn't get here.")
if PY3:
range_top = p
else:
# xrange on python2 can take integers representable as C long only
range_top = min(0x7fffffff, p)
for b in xrange(2, range_top):
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.")
if GMPY2:
def inverse_mod(a, m):
"""Inverse of a mod m."""
if a == 0:
return 0
return powmod(a, -1, m)
elif GMPY:
def inverse_mod(a, m):
"""Inverse of a mod m."""
# while libgmp likely does support inverses modulo, it is accessible
# only using the native `pow()` function, and `pow()` sanity checks
# the parameters before passing them on to underlying implementation
# on Python2
if a == 0:
return 0
a = mpz(a)
m = mpz(m)
lm, hm = mpz(1), mpz(0)
low, high = a % m, m
while low > 1:
r = high // low
lm, low, hm, high = hm - lm * r, high - low * r, lm, low
return lm % m
else:
def inverse_mod(a, m):
"""Inverse of a mod m."""
if a == 0:
return 0
lm, hm = 1, 0
low, high = a % m, m
while low > 1:
r = high // low
lm, low, hm, high = hm - lm * r, high - low * r, lm, low
return lm % m
try:
gcd2 = math.gcd
except AttributeError:
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, integer_types)
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