#! /usr/bin/python
# -*- coding: UTF-8 -*-
# formatted for tab-stops 4
#
# By Peter Schorn, peter.schorn@acm.org, December 2004 - November 2005
#
# Decompose a non-negative integer into the sum of at most four squares.
# Last change 25-Nov-2005, simplified decomposeProbablePrime
# Last change 13-Nov-2005, optimized search loop
# Last change 30-Oct-2005, improved isqrt function
# Last change 8-May-2005, improved isqrt function
import sys, operator, math, random
# Return the greatest common divisor using Euclid's algorithm
def gcd(a, b):
while b:
a, b = b, a % b
return a
# 'l' is a list of tuples where the first elements form the argument to 'func'
# and the last element is the expected result. Returns True iff all tests
# succeed.
def testfunc(func, l):
return False not in [func(*x[:-1]) == x[-1] for x in l]
# Test the gcd function
# Returns True if no error has been detected
def testgcd():
return testfunc(gcd, [(4, 6, 2), (6, 4, 2),
(1, 1, 1), (2, 1, 1), (1, 2, 1), (13000, 17000, 1000), (3, 4, 1)])
# Compute the jacobi symbol (a / b)
# Precondition: b odd
# Postcondition: Result is jacobi symbol (a / b) or 0 iff gcd(a, b) != 1
def jacobi(a, b):
s = 1
while a > 1:
d4 = a >> 2
m4 = a & 3
if m4: # a % 4 != 0
if m4 == 2: # a % 4 = 2 (i.e. a % 2 = 0)
if b & 7 in [3, 5]:
s = -s
a >>= 1
else:
if m4 == b & 3 == 3:
s = -s
a, b = b % a, a
else: # a % 4 = 0
a = d4 # a = a / 4
return a and s
# Test the jacobi function by comparing with a pre-computed table
# Returns True if no error has been detected
def testjacobi():
return testfunc(jacobi, [
(1, 1, 1), (1, 3, 1), (1, 5, 1), (1, 7, 1), (1, 9, 1), (1, 11, 1),
(1, 13, 1), (1, 15, 1), (1, 17, 1), (1, 19, 1), (1, 21, 1), (1, 23,
1), (1, 25, 1), (1, 27, 1), (1, 29, 1), (1, 31, 1), (2, 1, 1), (2, 3,
-1), (2, 5, -1), (2, 7, 1), (2, 9, 1), (2, 11, -1), (2, 13, -1), (2,
15, 1), (2, 17, 1), (2, 19, -1), (2, 21, -1), (2, 23, 1), (2, 25, 1),
(2, 27, -1), (2, 29, -1), (2, 31, 1), (3, 1, 1), (3, 3, 0), (3, 5,
-1), (3, 7, -1), (3, 9, 0), (3, 11, 1), (3, 13, 1), (3, 15, 0), (3,
17, -1), (3, 19, -1), (3, 21, 0), (3, 23, 1), (3, 25, 1), (3, 27, 0),
(3, 29, -1), (3, 31, -1), (4, 1, 1), (4, 3, 1), (4, 5, 1), (4, 7, 1),
(4, 9, 1), (4, 11, 1), (4, 13, 1), (4, 15, 1), (4, 17, 1), (4, 19, 1),
(4, 21, 1), (4, 23, 1), (4, 25, 1), (4, 27, 1), (4, 29, 1), (4, 31,
1), (5, 1, 1), (5, 3, -1), (5, 5, 0), (5, 7, -1), (5, 9, 1), (5, 11,
1), (5, 13, -1), (5, 15, 0), (5, 17, -1), (5, 19, 1), (5, 21, 1), (5,
23, -1), (5, 25, 0), (5, 27, -1), (5, 29, 1), (5, 31, 1), (6, 1, 1),
(6, 3, 0), (6, 5, 1), (6, 7, -1), (6, 9, 0), (6, 11, -1), (6, 13, -1),
(6, 15, 0), (6, 17, -1), (6, 19, 1), (6, 21, 0), (6, 23, 1), (6, 25,
1), (6, 27, 0), (6, 29, 1), (6, 31, -1), (7, 1, 1), (7, 3, 1), (7, 5,
-1), (7, 7, 0), (7, 9, 1), (7, 11, -1), (7, 13, -1), (7, 15, -1), (7,
17, -1), (7, 19, 1), (7, 21, 0), (7, 23, -1), (7, 25, 1), (7, 27, 1),
(7, 29, 1), (7, 31, 1), (8, 1, 1), (8, 3, -1), (8, 5, -1), (8, 7, 1),
(8, 9, 1), (8, 11, -1), (8, 13, -1), (8, 15, 1), (8, 17, 1), (8, 19,
-1), (8, 21, -1), (8, 23, 1), (8, 25, 1), (8, 27, -1), (8, 29, -1),
(8, 31, 1), (9, 1, 1), (9, 3, 0), (9, 5, 1), (9, 7, 1), (9, 9, 0), (9,
11, 1), (9, 13, 1), (9, 15, 0), (9, 17, 1), (9, 19, 1), (9, 21, 0),
(9, 23, 1), (9, 25, 1), (9, 27, 0), (9, 29, 1), (9, 31, 1), (10, 1,
1), (10, 3, 1), (10, 5, 0), (10, 7, -1), (10, 9, 1), (10, 11, -1),
(10, 13, 1), (10, 15, 0), (10, 17, -1), (10, 19, -1), (10, 21, -1),
(10, 23, -1), (10, 25, 0), (10, 27, 1), (10, 29, -1), (10, 31, 1),
(11, 1, 1), (11, 3, -1), (11, 5, 1), (11, 7, 1), (11, 9, 1), (11, 11,
0), (11, 13, -1), (11, 15, -1), (11, 17, -1), (11, 19, 1), (11, 21,
-1), (11, 23, -1), (11, 25, 1), (11, 27, -1), (11, 29, -1), (11, 31,
-1), (12, 1, 1), (12, 3, 0), (12, 5, -1), (12, 7, -1), (12, 9, 0),
(12, 11, 1), (12, 13, 1), (12, 15, 0), (12, 17, -1), (12, 19, -1),
(12, 21, 0), (12, 23, 1), (12, 25, 1), (12, 27, 0), (12, 29, -1), (12,
31, -1), (13, 1, 1), (13, 3, 1), (13, 5, -1), (13, 7, -1), (13, 9, 1),
(13, 11, -1), (13, 13, 0), (13, 15, -1), (13, 17, 1), (13, 19, -1),
(13, 21, -1), (13, 23, 1), (13, 25, 1), (13, 27, 1), (13, 29, 1), (13,
31, -1), (14, 1, 1), (14, 3, -1), (14, 5, 1), (14, 7, 0), (14, 9, 1),
(14, 11, 1), (14, 13, 1), (14, 15, -1), (14, 17, -1), (14, 19, -1),
(14, 21, 0), (14, 23, -1), (14, 25, 1), (14, 27, -1), (14, 29, -1),
(14, 31, 1), (15, 1, 1), (15, 3, 0), (15, 5, 0), (15, 7, 1), (15, 9,
0), (15, 11, 1), (15, 13, -1), (15, 15, 0), (15, 17, 1), (15, 19, -1),
(15, 21, 0), (15, 23, -1), (15, 25, 0), (15, 27, 0), (15, 29, -1),
(15, 31, -1)])
# Return True iff n is a prime. Algorithm uses trial division and is therefore
# only usable for 'small' n. This function is used to test the probable prime
# functions which work for 'large' n as well.
def isPrime(n):
if n <= 1:
return False
if (n <= 3) or (n == 5):
return True
if (n & 1 == 0) or (n % 3 == 0):
return False
t = 5
d = 2
while (t * t <= n) and (n % t):
t += d
d = 6 - d # skip trial divisors t with t % 3 = 0
return n % t != 0
# Low quality primality test
# Precondition: n is odd
# For n > 2: if result is False then n is composite
# if result is True then n maybe prime
# 341 (= 11 * 31) is smallest composite number where function returns True
def isProbablePrime2(n):
return (n > 1) and (pow(2, n - 1, n) == 1)
# Test the probable prime detector 'func' on all odd numbers between
# n and upperN. Return False iff there is a case where a prime number
# is erroneously flagged as composite.
def testisProbablePrimeFunc(func, n, upperN):
if n & 1 == 0:
n += 1 # Make sure that n is odd
while n <= upperN:
if func(n) < isPrime(n): # Error if composite was prime after all
return False
n += 2 # Next odd
return True
# The list of all odd primes less than 100
oddPrimes97 = filter(isProbablePrime2, range(100))
# The product of all odd primes less than 100
oddPrimeProduct97 = reduce(operator.mul, oddPrimes97)
# Good quality primality test
# Precondition: n is odd
# For n > 2: if result is False then n is composite
# if result if True then n maybe prime
# 42799 (= 127 * 337) is smallest composite number where function returns True
# Idea is to check the identity jacobi(2, p) = 2 ^ ((p - 1) / 2) % p
# which holds for primes (according to the Euler criterion)
# Note that isProbablePrime(1093 * 1093) = True
def isProbablePrime(n):
return (n > 1) and ((n < 561) or gcd(oddPrimeProduct97, n) == 1) and \
(pow(2, n >> 1, n) == ((n & 7 in [1, 7]) or (n - 1)))
# Use the fact that True == 1
# Return True iff n is 2 or a strong pseudo prime to base a
# Checking for strong pseudo primality is the core of the Miller-Rabin
# probabilistic prime check: if a number is a strong pseudo prime for
# sufficiently many bases, it is highly likely to be a prime. Note that this
# idea does not work with just checking a^(n-1) = 1 mod n as a Carmichael
# number c (such as 561) will satisfy this criterion for all a with gcd(a,c) = 1
def isStrongPseudoPrime(a, n):
if n == 2:
return True
n1 = n - 1 # i.e. -1 % n
r = n1 >> 1
k = 1
q = r >> 1
m = r & 1
while m == 0:
k += 1
r = q
q = r >> 1
m = r & 1
# Postcondition: (n = 2^k * r + 1) and (r % 2 = 1)
t = pow(a, r, n)
pt = n1 # pt will hold previous value of t
while (t != 1) and (t != n1) and k:
pt = t # save previous value
t = (t * t) % n # square modulo n
k -= 1
# We have a strong pseudo prime iff
# Case 1) We reached 1 and the previous value was -1 (mod n) or
# Case 2) We reached -1 and we have at least one squaring to go
return (t == 1) and (pt == n1) or (t == n1) and (k > 0)
# Note: We cannot replace (k > 0) with k as a boolean value is
# expected.
primeProduct3 = 5*7*13*23*29*31*37*53*61*97*127*149*151*157*137*229
# Very good quality primality test (but slower than isProbablePrime)
# For n > 0: if result is False then n is composite or 1
# if result if True then n maybe prime
# 220729 (= 103 * 2143) is smallest composite number where function returns True
def isProbablePrime3(n):
return (n > 1) and ((n < 2047) or gcd(primeProduct3, n) == 1) and \
isStrongPseudoPrime(2, n)
# Test the probable prime detector 'func' on all odd numbers between
# n and upperN. Return False iff there is a discrepancy between
# 'func' and isPrime
def testisProbablePrimeFuncStrict(func, n, upperN):
if n & 1 == 0:
n += 1 # Make sure that n is odd
while n <= upperN:
if func(n) != isPrime(n): # Error if functions disagree
print 'Prime error for %i. Got %s which is not correct.' \
% (n, func(n))
return False
n += 2 # Next odd
return True
# Probable prime check to be used. Note that correctness of algorithms in this
# module only depends on the fact that composite numbers are detected correctly.
useIsProbablePrime = isProbablePrime # productive use
#useIsProbablePrime = isProbablePrime # productive use
#useIsProbablePrime = isProbablePrime2 # testing
#useIsProbablePrime = lambda x: x > 1 # testing
#useIsProbablePrime = isProbablePrime3 # testing
# Performance evaluation
# Case 'small': time decompose.py 0 100000 f
# Case 'medium': time decompose.py 'sys.maxint-500L' 100000 f
# Case 'large': time decompose.py '10**200' 50 f
#
# (Mac OS X 10.3.8, 1.5GHz MHz PowerPC G4, 256 KB L2 Cache, 2MB L3 Cache,
# 1.5GB Memory, Python 2.3 (#1, Sep 13 2003, 00:49:11)
# [GCC 3.3 20030304 (Apple Computer, Inc. build 1495)] on darwin )
# small large total
# isProbablePrime 0:11.47 0:05.45 16.92
# isProbablePrime2 0:10.08 0:27.89 37.97
# lambda x: x > 1 0:10.90 0:27.73 38.63
# isProbablePrime3 0:13.44 0:13.66 27.10
# Cache for the odd primes less than 500
oddPrimesCache = filter(isProbablePrime, range(500))
# Return the factorial of n
def fac(n):
return reduce(operator.mul, range(1, n + 1), 1)
# Return the product of the first n primes
# pp(0) = 1, pp(1) = 2, pp(2) = 6, ...
def pp(n):
if n > 0:
return reduce(operator.mul, [p for p in oddPrimes(n - 1)], 1) << 1
return 1
def testpp():
return testfunc(pp, [(-10, 1), (-1, 1), (0, 1), (1, 2), (2, 6), (3, 30),
(9, 223092870)])
# Return the smallest probable prime p>=n with p % (m/gcd(m,r))=(r/gcd(m,r))
def dp(n, m, r):
g = gcd(m, r)
m /= g
p = m * (n / m) + r / g
if p < n:
p += m
while not isProbablePrime(p):
p += m
return p
def testdp():
return testfunc(dp, [(100, 4, 1, 101), (101, 4, 1, 101), (100, 6, 4, 101),
(101, 6, 4, 101)])
# Return next probable prime >= n
# Precondition: n is odd
def nextProbablePrime(n):
while not useIsProbablePrime(n):
n += 2
return n
# Return next probable prime >= n
def np(n):
if n <= 2:
return 2
return nextProbablePrime((n & 1) and n or (n + 1))
# For given even k and odd prime p, return True if n = k * p + 1 is prime.
# If True is returned, n is guaranteed to be prime but note that the result
# might be False even though n is prime.
def isnpk(k, p):
n = k * p + 1
if gcd(n, oddPrimeProduct97) > 1:
return n in oddPrimes97 # in case n is a prime less than 100
t = pow(2, k, n)
if t == 1:
return False
return pow(t, p, n) == 1
def testisnpk():
return testfunc(isnpk, [(4, 13, True), (2, 5003, True), (4, 5003, False),
(2, (1L << 31) - 1, False), (46, (1L << 31) - 1, True), (32, 3, True),
(20, 5, True), (22, 5, False)])
# For given k and odd prime p return K >= k such that n = K * p + 1 is prime.
# Note that n is then guaranteed to be prime but K is not necessarily the
# smallest value >= k to make this happen. Note also that Dirichlet's theorem
# guarantees the existence of such a prime (but this algorithm might fail
# for all such primes ...)
def npk(k, p):
if k & 1:
k += 1 # k is now even
while not isnpk(k, p):
k += 2
return k
# Generator for the odd primes
def oddPrimes(limit = -1):
index = 0
while (limit < 0) or (index < limit):
if index >= len(oddPrimesCache):
oddPrimesCache.append(nextProbablePrime(oddPrimesCache[-1] + 2))
yield oddPrimesCache[index]
index += 1
# Probabilistic check whether n is a square. If the result is (False, q) then n
# cannot be a square and (q / n) = -1. If the result is (True, 0), then n may
# or may not be a square.
# The idea is to compute the jacobi symbol (q / n) for 'certainty' many odd
# primes q. Since (q / n1 * n2) = (q / n1) * (q / n2), the jacobi symbol (q / n)
# is always 1 if n is a square. This also means that if we can find a q such
# that (q / n) = -1 then n cannot be a square.
# Heuristically each additional q filters out above half of the non-squares.
# The number n = 182919 is the smallest non-square for which
# isSquareProbabilistic(n, 15) returns True, 0
# The number n = 8042479 is the smallest non-square for which
# isSquareProbabilistic(n, 30) returns True, 0
def isSquareProbabilistic(n, certainty):
generator = oddPrimes()
while certainty:
q = generator.next()
if jacobi(n % q, q) == -1: # q is a witness that n is not a square
return False, q
certainty -= 1
return True, 0
# There are exactly 336 squares mod 10080. This is the smallest rate of squares
# for all n < 15840.
magicN = 10080
squaresModMagicN = {}.fromkeys([(i * i) % magicN for i in range(magicN >> 1)])
# If result is True then n may or may not be a square. If the result is False
# then n cannot be a square. The 'false positive' rate is 3.3% = 336 / 10080 =
# 1 / 30. The smallest non-square for which isProbableSquare(n) = True is 385.
def isProbableSquare(n):
return (n % magicN) in squaresModMagicN
def testisProbableSquare(low, high):
while low <= high:
if not isProbableSquare(low) and (isqrt(low)[1] == 0):
return False
low += 1
return True
# Compute an imaginary unit modulo p
# Precondition: p % 4 = 1
# Postcondition: p is prime implies: Result is (x, False) and x^2 % p = -1
# The result can also be (x, True) and then x^2 = p
# Note that this algorithm might succeed even when the argument is not prime.
# Example (1)
# n = 3277 = 29 * 113 % 8 = 5, 2^((n-1)/4) = 2^819 = 128 mod n
# 128^2 % n = -1
# Explanation:
# order(2) mod 29 = order(2) mod 113 = order(2) mod n = 28 divides n - 1
# This implies: 2^(n-1) mod 29 = 2^(n-1) mod 113 = 1
# We also have 2^((n-1)/2) mod 29 = 2^((n-1)/2) mod 113 = -1
# overall 2^((n-1)/2) % n = -1
# Example (2)
# p = 3281 = 17 * 193 % 8 = 1, jacobi(3, p) = -1
# 3^((p-1)/4) = 3^820 = 81 mod p, 81^2 = -1 mod p
# Example (3)
# p = 49141 = 157 * 313 % 8 = 5, 2^((p-1)/4) = 2^12285 = 32527 mod p
# 32527^2 = -1 mod p, isProbablePrime(p) = True
# Example (4)
# p = 665281 = 577 * 1153 % 8 = 1, jacobi(13, p) = -1
# 13 ^ ((p-1)/4) = 13 ^ 166320 = 531393 mod p
# 531393^2 = -1 mod p, isProbablePrime(p) = True
# Note also that if p is a square, jacobi(q, p) = 1 for all q
def iunit(p):
if p & 7 == 5:
q = 2
else:
primes = oddPrimes()
q = primes.next() # q is an odd (probable) prime
# (q % 2 = 1) and (p % 4 = 1) implies jacobi(q, p) = jacobi(p % q, q)
while jacobi(p % q, q) == 1: # jacobi(q, p)
q = primes.next()
if (q == 229) and isProbableSquare(p):# loop is running quite long
# reached for p = dp(1, 4*pp(k), 1) for k > 47 or p = 1093 ** 2
s, r = isqrt(p) # check, if p is a square
if r == 0: # if yes, return square root
return s, True
return pow(q, p >> 2, p), False
# Test the iunit function for arguments between p and upperP making sure
# that the argument is a (probable) prime congruent 1 modulo 4
# Returns True if no error has been detected
def testiunit(p, upperP):
p = ((p >> 2) << 2) + 1 # Make sure that p % 4 = 1
upperP = min(upperP, 42799) # and that upperP does not exceed range where
while p <= upperP: # isProbablePrime works correctly.
if isProbablePrime(p):
if (iunit(p)[0] ** 2 + 1) % p: # This is the check
return False
p += 4 # maintain invariant p % 4 = 1
return True
ln2 = math.log(2)
# Returns the smallest integer greater or equal to the logarithm of n to base 2,
# i.e. the result r satisfies 2^r <= n < 2^(r+1).
def floorLog2(n):
r = int(0.5 + math.log(n) / ln2)
return ((1L << r) <= n) and r or (r - 1)
# The square root algorithm uses its iterative variant for numbers less than
# 2^1024. For larger numbers the recursive case is invoked.
isqrtIterativeBetter = 1L << 1024
# Input: n >= 0, log2n = floor(log(2, n)), 2^log2n <= n < 2^(log2n + 1)
# Output: (s, r) where (s-1)^2 <= n < (s+1)^2 and s^2 + r = n
# For very small n, the built-in square root function is invoked. For larger n
# (typically n < 10^300) an iterative algorithm is used while for the largest
# n, a divide & conquer approach is employed almost identical to the one from
# Paul Zimmermann (Karatsuba Square Root, http://www.inria.fr/rrrt/rr-3805.html)
# The difference is that a slightly weaker invariant ((s-1)^2 <= n < (s+1)^2
# can be shown to suffice - note that the r = n - s^2 can be negative in the
# case that the square root is over estimated.
def isqrtInternal(n, log2n):
if n <= sys.maxint:
s = int(math.sqrt(n))
return s, n - s * s
if n < isqrtIterativeBetter:
d = 7 * (log2n / 14 - 1)
q = 7
# assert (0 < (n >> (d + d)) < (1 << 28)) and (q <= d)
s = int(math.sqrt(n >> (d + d)))
while d:
# assert (s - 1) ** 2 <= (n >> (d + d)) < (s + 1) ** 2
if q > d:
q = d
s <<= q
d -= q
q += q
s = (s + (n >> (d + d)) / s) >> 1
return s, n - s * s
log2b = log2n >> 2
mask = (1L << log2b) - 1L
s, r = isqrtInternal(n >> (log2b + log2b), log2n - (log2b + log2b))
q, u = divmod((r << log2b) + ((n >> log2b) & mask), s + s)
return (s << log2b) + q, (u << log2b) + (n & mask) - q * q
# Input: n >= 0
# Output: (s, r) where s^2 <= n < (s+1)^2 and s^2 + r = n
def isqrt(n):
assert n >= 0
if n <= sys.maxint:
s = int(math.sqrt(n))
return s, n - s * s
s, r = isqrtInternal(n, floorLog2(n))
return (r < 0) and (s - 1, r + s + s - 1) or (s, r)
# Test the isqrt function for arguments between n and upperN
# Returns True if no error has been detected
def testisqrt(n, upperN):
while n <= upperN:
s, r = isqrt(n)
if not ((s * s <= n < (s + 1) ** 2) and (s * s + r == n)):
return False # return False iff postcondition of isqrt not met
n += 1
return True
# Try to decompose n into a sum of two squares
# Precondition: n is a probable prime, n % 4 = 1
# (a, b, True) with a^2 + b^2 = n iff decomposition was successful
# (0, 0, False) iff decomposition was not successful
def decomposeProbablePrime(n):
b, r = iunit(n)
if r:
return 0, b, True
if (b * b + 1) % n: # Check whether we got an imaginary unit
return 0, 0, False # Indicate failure if not
a = n
while b * b > n:
a, b = b, a % b
return b, a % b, True
# Perform n tests starting at the smallest prime p >= start with p % 4 = 1
def testdecomposeProbablePrime(start, n):
p = ((start >> 2) << 2) + 1
if p < start:
p += 4
while n:
while not isProbablePrime(p):
p += 4
a, b, c = decomposeProbablePrime(p)
if c and (a * a + b * b != p):
return False
n -= 1
p += 4
return True
# Dictionary of exceptional cases the search in decompose cannot handle
decomposeExceptions = {
9634: [56, 57, 57], 2986: [21, 32, 39], 1906: [13, 21, 36],
1414: [ 6, 17, 33], 730: [ 0, 1, 27], 706: [15, 15, 16],
526: [ 6, 7, 21], 370: [ 8, 9, 15], 226: [ 8, 9, 9],
214: [ 3, 6, 13], 130: [ 0, 3, 11], 85: [ 0, 6, 7],
58: [ 0, 3, 7], 34: [ 3, 3, 4], 10: [ 0, 1, 3],
3: [ 1, 1, 1], 2: [ 0, 1, 1] }
# Decompose an integer n>=0 into a sum of up to four squares. Result is a
# tuple where the first entry is a list of four integers and the second
# entry is True iff the search loop was successful. Note that currently
# there is no known input where the search loop would not be successful
# but there is no known proof.
def decompose(n):
# Multiply all elements of l with v and return a tuple with the sorted list
# and the value True indicating success
def sortScaled(l):
l.sort()
return [v * x for x in l], True
assert n >= 0
if n <= 1: # Done if n = 0 or n = 1
return [0, 0, 0, n], True
v = 1 # Remove powers of 4
while n & 3 == 0: # n % 4 = 0
n >>= 2 # n = n / 4
v += v
# Postcondition: (nOld = v^2*n) and (n%4 != 0) and (v = 2^k for some k >= 0)
sq, p = isqrt(n)
if p == 0: # n is a square
return [0, 0, 0, v * sq], True
if (n & 3 == 1) and useIsProbablePrime(n): # n is prime, n % 4 = 1
s, r, done = decomposeProbablePrime(n)
if done: # otherwise n was not a prime
return sortScaled([0, 0, s, r])
if n & 7 == 7: # Need 4 squares, subtract largest square delta^2
# such that (n > delta^2) and (delta^2 % 8 != 0)
delta, n = sq, p
if sq & 3 == 0:
delta -= 1
n += delta + delta + 1
# sq % 4 = 0 -> delta % 8 in [3, 7] ->
# delta^2 % 8 = 1 -> n % 8 = 6
# sq % 4 != 0 -> delta % 8 in [1, 2, 3, 5, 6, 7] ->
# delta^2 % 8 in [1, 4] -> n % 8 in [3, 6]
# and this implies n % 4 != 1, i.e. n cannot be a sum of two squares
sq, p = isqrt(n)
# sq^2 != n since n % 4 = 3 which is never true for a square
else:
delta = 0
# Postcondition: (sq = isqrt(n)) and (n % 8 != 7) and (n % 4 != 0) and
# (nOld = v^2 * (n + delta^2))
# This implies that n is a sum of three squares - now check whether n
# is one of the special cases the rest of the algorithm could not handle.
if n in decomposeExceptions:
# Retrieve pre-computed result and scale with v
return sortScaled([delta] + decomposeExceptions[n])
# Now perform search distinguishing two cases noting that n % 4 != 0
# Case 1: n % 4 = 3, n = x^2 + 2*p, p % 4 = 1,
# p prime, p = y^2 + z^2 implies n = x^2 + (y + z)^2 + (y - z)^2
if n & 3 == 3:
if sq & 1 == 0:
sq -= 1
p += sq + sq + 1
p >>= 1
while True:
if useIsProbablePrime(p):
r0, r1, done = decomposeProbablePrime(p)
if done:
return sortScaled([delta, sq, r0 + r1, abs(r0 - r1)])
sq -= 2 # Next sq
if sq < 0: # Should not happen, no known case
return 4 * [0], False
p += sq + sq + 2 # Next p
# Case 2: n % 4 = 1 or n % 4 = 2, n = x^2 + p,
# p % 4 = 1, p prime, p = y^2 + z^2 implies n = x^2 + y^2 + z^2
if (n - sq) & 1 == 0:
sq -= 1
p += sq + sq + 1
while True:
if useIsProbablePrime(p):
r0, r1, done = decomposeProbablePrime(p)
if done:
return sortScaled([delta, sq, r0, r1])
sq -= 2 # Next sq
if sq < 0: # Should not happen, no known case
return 4 * [0], False
p += 4 * (sq + 1) # Next p
# Decompose all integers from list li into at most four squares and print
# the results iff verbose = True. The result is a list of all cases where the
# decomposition was not successful. Note that this list should be empty.
def checkDecompose(li, verbose):
f = []
for i in li:
r, done = decompose(i)
if verbose: # Use the fact that r is sorted
if r[0]: # Need four squares
s = '%i = %i^2 + %i^2 + %i^2 + %i^2' % (i, r[0], r[1],
r[2], r[3])
elif r[1]: # Need three squares
s = '%i = %i^2 + %i^2 + %i^2' % (i, r[1], r[2], r[3])
elif r[2]: # Need two squares
s = '%i = %i^2 + %i^2' % (i, r[2], r[3])
else: # Number is a square
s = '%i = %i^2' % (i, r[3])
# print also number of squares needed
print '%s # %i' % (s, 4 - r[:-1].count(0))
if not done or (sum([x*x for x in r]) != i):
f.append((i, r, done))
return f
# Perform a small self test of the major function.
# Return True iff all tests successful.
def selfTest():
return (testgcd() and
testisProbablePrimeFunc(isProbablePrime2, 1, 1000) and
testisProbablePrimeFunc(isProbablePrime, 1, 1000) and
testjacobi() and
testisqrt(0, 125) and
testisqrt(10**20 + random.randint(0, 10**19), 20) and
testiunit(5, 2377) and
(iunit(1093 ** 2)[0] == 1093) and
testisProbablePrimeFuncStrict(isProbablePrime3, 3, 1000) and
testisProbablePrimeFuncStrict(isProbablePrime, 3, 1000) and
testisnpk() and
testpp() and
testdecomposeProbablePrime(5, 100) and
testdecomposeProbablePrime(10 ** 5 + 1, 3) and
testisProbableSquare(100, 400) and
(len(checkDecompose([dp(1, 4 * pp(48), 1)], False)) == 0) and
testdp()
)
def main():
# Print usage information for the script and exit.
def usage():
print '''\
Usage: %s start(Exp) [iterations(Exp)] [v(char)] [f(char)]
Decomposes the integers in start (if start evaluates to a list) or in
[start, start + iterations - 1] (if start evaluates to an integer)
into the sum of at most 4 squares.
Expression can contain the following functions:
gcd(a,b): greatest common divisor of a and b
jacobi(a,b): the Jacobi symbol of a and b
fac(n): factorial of n
pp(n): product of the first n primes [pp(0)=1, pp(1)=2, pp(2)=6, ...]
dp(n,m,r): smallest probable prime p >= n with p mod m/gcd(m,r)=r/gcd(m,r)
np(n): smallest probable prime p >= n
isqrt(n)[0]: integer square root of n is y with y^2 <= n < (y + 1)^2
v: verbose f: fast (omit self test)
Examples:
decompose.py 10007
decompose.py '577*1153'
decompose.py 'np(19081961*10**3)' 2 v
decompose.py '[2*10**i+19081961 for i in [34, 35, 36]]' v
time decompose.py 0 100000
''' % sys.argv[0]
sys.exit(2)
# Called because some form of command line syntax error detected
# Return the value of 'intString' if it contains a valid expression.
# Otherwise print 'errorMessage', the usage information and exit.
def getInt(intString, errorMessage):
try:
return eval(intString)
except:
print errorMessage % intString
usage()
# Return True iff string 'opt' is in the lower case version
# of the script's arguments.
def hasOption(opt):
return opt in [x.lower() for x in sys.argv]
if not (2 <= len(sys.argv) <= 5):
usage()
if not (hasOption('f') or selfTest()):
print 'Self test failed! Exiting...'
sys.exit(1)
li = getInt(sys.argv[1], 'Illegal expression \'%s\' for start.')
verbose = hasOption('v')
if not isinstance(li, list):
# Create list of numbers to check.
# If number of iterations is 0, use 1 instead.
li = range(li, li + ((len(sys.argv) > 2) and
not sys.argv[2].lower() in ['f', 'v'] and
getInt(sys.argv[2],
'Illegal expression \'%s\' for number of iterations.') or 1))
r = checkDecompose(li, verbose or (len(li) == 1))
if len(r):
print 'Errors: ', r
if __name__ == '__main__':
main()