#
# New things are on the top...
#
# Paul Garrett, 1998-...
#
import math # use sqrt, log, exp, ceil
import __main__
######################################################################
def lg( n ):
return math.log(0.0+n)/math.log(0.0+2)
def hamming_weight( n ):
wt = 0
while n > 0:
wt = wt + ( n % 2 )
n = n/2
return wt
def wt( n ):
return hamming_weight( n )
def entropy( p ):
return -p * lg(p) - (1-p)*lg(1-p)
def H(p):
return entropy(p)
def channel_capacity( p ):
return 1 - H(p)
def cap(p):
return channel_capacity(p)
#####################################################################
def is_prime( n ):
if (factors(n) == [n]):
return 1
else:
return 0
def binomial_coefficient( n , k ):
out = n + 0L
for i in range(1,k):
out = out * (n-i)
for i in range(1,k+1):
out = out/i
return out
def inverse( x, n ):
""" Usage: inverse( x, n )
computes inverse of x modulo n, assuming gcd(b,n)=1
via Euclidean algorithm
If n=0, returns rational number inverse..."""
if n != 0:
n = abs(n)
return (full_euclid(n,x)[1] + n) % n
else:
return 1.0/x
def sgn( n ):
if n > 0:
return 1
elif n < 0:
return -1
else:
return 0
def full_euclid( x,y):
#
# Hmmm. Note that divmod(17,-3) == (-6,-1) Rrrr...
#
""" Usage: full_euclid( x, y )
Returns [a,b,g] so that a.x + b.n = g = gcd(x,n) """
a,b,c,d = 1L,0L,0L,1L
sgnx = sgn(x)
sgny = sgn(y)
x = abs(x)
y = abs(y)
while y != 0:
[q,r] = divmod(x,y)
x,y = y,r
a,b,c,d = c,d,a-c*q,b-q*d
# now a.x_abs + b.y_abs = gcd
return a*sgnx, b*sgny, x # now x is the gcd
def sqrt( n ):
""" Usage: sqrt( n )
Is infinite-precision long-integer square-root evaluator
Uses Newton's method, not floating-point """
if n < 1000000:
return int(math.sqrt(n))
digits = len( repr(n) ) - 1
if digits <= 0: digits = 1
x = 10L ** (digits/2) ## smaller than actual sqrt
decr = (x - n/x)/2
while decr >= 1 or decr <= -1:
x = x - decr
decr = (x - n/x)/2
return x
######################################################################
def find_last_index( bound, List ):
""" Usage: find_index_in_ordered_list( bound, List )
finds highest index i so that List[i] <= bound
Uses divide-and-conquer (binary expansion) search """
lg = int( math.log( len(List)-1 )/ math.log( 2 ) )
now = pow(2, lg) - 1
while ( lg > 0 and (List[now] > bound or
( now+1 < len(List) and List[now+1] <= bound ))):
lg = lg - 1
inc = pow(2,lg)
if List[now] > bound: # guess is too high
now = now - inc
else: # too low
if now + inc < len(List):
now = now + inc
else:
pass
return now
######################################################################
def _mr(power_of_2, m, b, n):
# private function called by "public" miller_rabin()
t = 0
m = m + 0L
n = n + 0L
b = b + 0L
Y = pow( b, m, n)
if (Y == 1):
return 1
else:
while ( t < power_of_2 ):
if (Y == n-1):
return 1
Y = Y * Y % n
if (Y == 1):
return 0
t = t + 1
return 0
def standard( n ): ## commercial pseudoprimality test
bases_list = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,73]
return miller_rabin( n , bases = bases_list )
#
# In the latter, there seems to be no point in trial division first.... !?!?
#
# (as indicated by naive benchmarking tests on 100-digit numbers)
#
def miller_rabin( n , bases = [2] ):
""" Usage: miller_rabin( n, bases = [2] )
Performs Miller-Rabin test on n with list "bases"
returns 0 for composite, 1 for probable prime """
if n < 2 :
n = 2
m = n-1L
n = n + 0L
power_of_2 = 0L
while (m % 2 == 0):
power_of_2 = 1 + power_of_2
m = m/2
for b in bases:
if ( _mr(power_of_2, m, b, n) == 0 ):
return 0
return 1
def factors( m, verbose = 0 ):
""" Usage: factors( m, verbose = 0 )
returns [list_of_factors]"""
list_of_factors = []
sq = sqrt(m)
while m % 2 == 0:
list_of_factors.append( 2 )
m = m/2
sq = sqrt(m)
if verbose:
print 2
d = 3
while d <= sq:
if m % d == 0:
list_of_factors.append( d )
m = m/d
sq = sqrt(m)
if verbose:
print d
else:
d = d+2
if m > 1:
list_of_factors.append( m )
if verbose:
print m
return list_of_factors
def gcd(x,y):
""" Usage: gcd(x,y) \n returns gcd"""
while ( y>0 ):
x = x % y
(x,y) = (y,x)
return x
def trial_division( n, bound = 1000):
""" Usage trial_division( bound, n )
returns smallest divisor d with 1 < d <= n, or 1 if none such"""
bound = bound + 2L
if ( n % 2 == 0 ):
return 2
d = 3L
while ( d <= bound ):
if ( n % d == 0 ):
return d
d = d + 2
return 1
def trial_division_by_list( n, the_list ):
""" Usage: trial_division_by_list( n, the_list )
returns smallest divisor >1 on the list, or 1 """
for tester in the_list:
if n % tester == 0:
return tester
return 1
def naive( n ):
""" Usage: naive( n )
returns smallest factor <= sqrt(n), or 1 if none such"""
sq = sqrt( n ) + 0L
return trial_division( n, sq )
def sophie_germain( n, ptest = naive ): # "ptest" should be a primality test
""" Usage: sophie_germain( n, ptest = naive )
finds next sophie germain prime above n (so that 2n+1 also prime)
using given primality test, with default "naive" test"""
if n % 2 == 0:
n = n + 1
while ( ptest( n ) == 0 or ptest( 2*n+1 ) == 0 ):
n = n+2
return [n, 2*n+1]
def fermat( n, bases = [2] ):
""" Usage: fermat( n, bases = [2] ) with n odd
returns 0 for composite, 1 for probably prime"""
sq = sqrt ( n ) ######## this is local, long-valued "sqrt"
for b in bases:
if ( b > sq ):
return 1
b = 0L + b
if gcd(b, n) > 1:
pass ################# else "carmichael" becomes vacuous!
elif pow(b, n, n) != b:
return 0
return 1
def _rho_random_function(x,n):
""" Usage: _rho_random_function(x,n)
it is x*x+1 % n """
return (x * x + 1L) % n
def rho(n, bound = 0, verbose = 0 ):
""" Usage: rho(n, bound = 0, verbose = 0 )
This is Pollard's rho method with some specific choices:
uses _rho_random_function() which by default is x^2+1 % n
prints number of cycles every 10,000 if verbose == true
if bound is set >0 then limits number of cycles before return
Returns [factor, cycle-count]"""
y = 2L ## edit value here?!
x = _rho_random_function( y,n )
count = 0
if verbose == 0 and bound == 0:
while 1:
count = count + 1
g = gcd(n, x-y)
if g > 1 and g < n:
return [g, count]
else: # may never return...
x = _rho_random_function( _rho_random_function( x,n ), n)
y = _rho_random_function( y,n )
elif verbose == 0 and bound > 0:
while count < bound:
count = count + 1
g = gcd(n, x-y)
if g > 1 and g < n:
return [g, count]
else: # may never return...
x = _rho_random_function( _rho_random_function( x,n ), n)
y = _rho_random_function( y,n )
return [1,count]
elif verbose == 1 and bound == 0:
while 1:
count = count + 1
g = gcd(n, x-y)
if g > 1 and g < n:
print " count = " + `count`
print g
return [g, count]
else:
x = _rho_random_function( _rho_random_function( x,n ), n)
y = _rho_random_function( y,n )
if (count % 10000 == 0): # this block is missing in non-verbose
print "count = " + `count`
count = count + 1
else: # verbose and bounded case
while count < bound:
count = count + 1
g = gcd(n, x-y)
if g > 1 and g < n:
print " count = " + `count`
print g
return [g, count]
else:
x = _rho_random_function( _rho_random_function( x,n ), n)
y = _rho_random_function( y,n )
if (count % 10000 == 0): # this block is missing in non-verbose
print " count = " + `count`
count = count + 1
print " No factor found in " + `bound` + " tries "
return [1, count]
class small_primes:
""" Usage: small_primes( bound )
Initializes object "sp" so that len(sp) = number primes under bound
and sp[i] retrieves the i-th prime (under the bound) """
def __init__(self, bound):
self.plist = []
self.plist.append( 2 )
t = 3
sq = sqrt(bound)
while (t < bound):
if is_prime(t):
self.plist.append( t )
t = t + 2
def get( self ):
return self.plist
def __len__( self ):
return len( self.plist )
def __getitem__( self, i ):
return self.plist[i]
######################################################################
def carmichael(lower=500L, upper=3000L, verbose = 0):
List = []
if lower < 100:
lower = 500L
primes_bound = sqrt( upper )
sp = small_primes( primes_bound )
# primes_bound = long(math.ceil( math.exp( math.log(upper)/3 ) ) )
# th = small_primes( primes_bound )
#
# Note: since a carmichael number is divisible by at least 3 distinct
# primes (!), we need only search for prime factors up through the
# cube root, rather than square root.
#
# This would have ramifications for rho-test, too, I guess... ?
#
if ( lower % 2 == 0 ):
lower = lower + 1L
else:
lower = lower + 0L
for t in xrange(lower, upper, 2L):
if ( fermat( t, [2L]) == 1 and
is_prime( t ) == 0 and
fermat( t, sp.plist ) == 1):
List.append(t)
if verbose:
print t
else:
pass
print List
######################################################################
#
# Lucas-Lehmer test for Mersenne numbers to be prime
#
def lucas_lehmer( s ):
""" Usage lucas_lehmer( s )
Applies Lucas-Lehmer test to assess primality of 2^s-1"""
if s == 1: return 0
if s == 2: return 1
n = pow(2L, s) - 1
for t in xrange(2, sqrt(s)+1, 2):
if s % t == 0: return 0 # composite exponent
u = 4
for i in xrange(0,s-2):
u = (u * u - 2) % n
if u == 0:
return 1
else:
return 0
def man( obj ):
""" Usage: man( obj )
Prints the document string obj.__doc__ attached to an object,
function, method, ... named "obj" """
print obj.__doc__
def is_primitive_root( modulus, b = 2L ):
""" Usage: is_primitive_root( modulus, b = 2 )
Returns 1 if b is a primitive root modulo, else 0"""
less_one = modulus - 1
list_of_factors = distinct_factors( modulus - 1 )
for factor in list_of_factors:
if pow(b, less_one/factor, modulus) == 1:
return 0
return 1
def primitive_root( modulus ):
""" Usage: primitive_root( modulus )
Returns smallest primitive root modulo modulus"""
b = 2
while ( is_primitive_root( modulus, b) == 0):
b = b+1
return b
def next_prime( start, test = naive, increment = 2 ):
""" Usage: next_prime( start, test = naive )
Returns next prime above "start"
Uses primality test specified, defaults to "naive",
which is simply trial division
"""
if ( start % 2 == 0 ):
start = start + 1
else:
start = start + increment
test_outcome = test( start )
while ( test_outcome == 0 or (test_outcome < start and test_outcome > 1 )):
start = start + increment
test_outcome = test( start )
return start
def p_minus( n, bound = 100000 ):
""" Usage p_minus( n, bound = 100000 )
Applies Pollard's p-1 test with smoothness bound "bound"
Returns factor, which is proper if successful """
plist = [2,3,5] ## ,5,7,11,13,17,19,23,29,31]
# sp = small_primes( bound+1 )
a = 3
print ("n = " + `n`)
print ("Use a = " + `a`)
for p in plist: ## sp.plist:
expr = long( math.log( n )/math.log( p ) )
print("exp for " + `p` + " is " + `expr`)
a = pow( a, pow(p+0L, expr), n)
print("a^that mod " + `n` + " is " + `a`)
#
# for t in xrange(0, expr+1): #### don't eval gcd's so often!?
# a = pow( a, p, n)
#
g = gcd( n, a - 1 )
print("gcd(n=" + `n` + ", a-1=" + `a-1` + ") is " + `g`)
if (g > 1 and g < n):
return g
return gcd( n, a - 1 )
def small_factors( n, bound = 1000000, verbose = 0 ): ## 29 Dec 98
""" Usage: small_factors( n, bound = 1000000 )
If stdout is unbuffered (with -u switch in interpreter) then
factors under bound are printed as they're found.
Return value is compound list
[ list-of-factors, leftover ] """
sq = sqrt( n )
if sq < bound:
bound = sq
List = []
for p in __main__.primes:
if p > bound or n == 1:
if verbose:
print "\n got all the factors",
break
while n % p == 0:
n = n/p
List.append(p)
sq = sqrt( n )
if sq < bound:
bound = sq
if verbose:
print p,
return [List, n]
def small_attack( n, cycle_bound = 200000, verbose = 1 ):
""" Usage: small_attack( n, cycle_bound = 200000, verbose = 1 )
If stdout is unbuffered (with -u switch in interpreter) then
factors are printed as they're found.
Return value is compound list
[ list-of-factors, leftover, 0 or 1 ]
where last indicates probable prime or not
First uses list of primes under 1,000,000, then
Pollard's rho test for "cycle_bound" cycles """
if verbose:
print "... doing trial division ..."
[smallFactors, n] = small_factors( n, 1000000, verbose )
if n == 1:
return [smallFactors, n, 1]
if miller_rabin( n, [2,3,5,7,11,13,17,19,23]) == 1:
if verbose:
print "\n ... didn't do rho ... "
return [smallFactors, n, 1]
if verbose:
print "\n ... doing rho now ..."
[div, cycles_used] = rho( n, cycle_bound, verbose )
n = n/div
if verbose == 0: # non-verbose case
while div > 1 and miller_rabin( n, [2,3,5,7,11,13,17,19,23]) == 0:
smallFactors.append( div )
[div, cycles_used] = rho( n, cycle_bound)
n = n / div
probable_prime = miller_rabin( n, [2,3,5,7,11,13,17,19,23] )
return [smallFactors, n, probable_prime]
else: # verbose case
while div > 1 and miller_rabin( n, [2,3,5,7,11,13,17,19,23]) == 0:
smallFactors.append( div )
[div, cycles_used] = rho( n, cycle_bound, verbose )
n = n / div
probable_prime = miller_rabin( n, [2,3,5,7,11,13,17,19,23] )
return [smallFactors, n, probable_prime]
def distinct_factors( m, verbose = 0):
""" Usage: distinct_factors( m, verbose = 0 )
returns [list_of_factors]"""
list_of_factors = []
sq = sqrt(m)
while m % 2 == 0:
list_of_factors.append( 2 )
m = m/2
sq = sqrt(m)
if verbose:
print 2
d = 3
while d <= sq:
if m % d == 0:
list_of_factors.append( d )
while m % d == 0:
m = m/d
sq = sqrt(m)
if verbose:
print d
else:
d = d+2
if m > 1:
list_of_factors.append( m )
if verbose:
print m
return list_of_factors
######################################################################