There are three files, separated by ----- lines.
The first file defines a class for computations on rationals
The second file gives a very good random generator
The third file gives two classes for modular computations, including the
chinese remainder theorem.
--------------------------------------------------------------
rational.py
# routines for rational numbers
# written by Jurjen NE Bos
# defines class "rat" that does rational computations
# creation of a rational with rat(number) or rat(numerator, denominator)
# number may be int, long or even float
# e.g. rat(1.0/7.0)
# methods: all normal operations are allowed
# num() and den() produce numerator and denominator respectively
import rational
from math import floor
def gcd(a, b):
while a<>0:
a, b = b%a, a
return b
def simplify(n, d):
g = gcd(n, d)
n, d = n/g, d/g
if d<0: n, d = -n, -d
return n, d
precision = 64-4 # precision of the approximation in bits
def float2fract(x):
sign, value = cmp(x,0), abs(x)
margin = pow(2.0,-precision)
if value<margin: return sign, long(1/value)
if value>1/margin: return long(x), 1
margin = margin * value
n = value
t1,n1,t0,n0 = 1,0,long(floor(n)),1
while abs(float(t0)/float(n0)-value)>margin:
n = 1/(n-floor(n))
t1,n1,t0,n0 = t0, n0, t1+t0*long(floor(n)), n1+n0*long(floor(n))
return t0*sign, n0
class rat:
def __init__(self, *a):
if len(a)==1:
if type(a[0])==int_type or type(a[0])==long_type:
self.n, self.d = a[0], 1
elif type(a[0])==rat_type:
self.n, self.d = a[0].n, a[0].d
elif type(a[0])==float_type:
self.n, self.d = float2fract(a[0]) #
compute approximation
else: raise TypeError, 'only numbers can be
converted to rat'
else:
self.n, self.d = simplify(a[0], a[1])
def num(self): return self.n
def den(self): return self.d
def __repr__(self): return 'rat(' + `self.n` + ', ' + `self.d` + ')'
def __cmp__(self, other):
if self.n==other.n and self.d==other.d: return 0
if self.n*other.d>self.d*other.n: return 1
else: return -1
def __add__(self, other):
#we have to coerce ourselves
if type(other)<>rat_type:
if type(other)==int_type or type(other)==long_type:
return self+rat(other)
if type(other)==float_type: return float(self)+other
return rat(self.n*other.d+self.d*other.n, self.d*other.d)
def __sub__(self, other):
return rat(self.n*other.d-self.d*other.n, self.d*other.d)
def __mul__(self, other):
#we have to coerce ourselves
if type(other)<>rat_type:
if type(other)==int_type or type(other)==long_type:
return self*rat(other)
if type(other)==float_type: return float(self)*other
return rat(self.n*other.n, self.d*other.d)
def __div__(self, other): return rat(self.n*other.d, self.d*other.n)
def __pow__(self, other):
exp=int(other)
if exp>=0: return rat(pow(self.n, exp), pow(self.d, exp))
return rat(pow(self.d,-exp), pow(self.n,-exp))
def __neg__(self): return rat(-self.n, self.d)
def __pos__(self): return self
def __abs__(self): return rat(abs(self.n), self.d)
def __nonzero__(self): return self.n<>0
def __int__(self):
if self.d<>1: raise ValueError, 'rat is not a whole number'
return int(self.n)
def __long__(self):
if self.d<>1: raise ValueError, 'rat is not a whole number'
return long(self.n)
def __float__(self): return(float(self.n)/float(self.d))
def __coerce__(self, other):
if type(other)==int_type or type(other)==long_type: return
self, rat(other)
if type(other)==float_type: return self.__float__(), other
def __hash__(self):
if self.d==1: return hash(n)
return hash(d)^hash(n)
int_type=type(0)
long_type=type(0L)
rat_type=type(rat(0))
float_type=type(0.0)
--------------------------------------------------------------
random.py
# high-quality random routine
# written by Jurjen N.E. Bos (jurjen@q2c.nl)
# passes the spectral test
# multiplier by C.E.Haynes (see Knuth, vol.2)
import random
from time import time
seed = 0x4A75726A656EL
modulus = 1L<<64
modulusMask = modulus-1
multiplier = 6364136223846793005L
def rdz(*s): # randomize from number. rdz() uses current time
if s: random.seed = long(s[0])
else: random.seed = long(time())
def next(): # generate next element from random sequence.
Internal use only
random.seed = (seed*multiplier+1)&modulusMask
return seed
# auxiliary numbers for rand
modulusSize = 64 # bits of modulus
resolution = 16 # number of reliable bits of random
complement = 48 # remaining bits from random number
def rand(*n): # rand(n): random integer 0..n-1; rand(): random
float 0..<1
# works OK if n very large
if n:
n = n[0]
# create a fraction with enough accuracy
num, den = 0L, 1L<<resolution
while den<n:
num = (num<<resolution)+(next()>>complement)
den = den<<resolution
num = (num<<modulusSize)+next()
den = den<<complement
return n*num/den
else:
return float(next())/float(modulus)
def choice(s): return s[int(rand(len(s)))] # random element from sequence
--------------------------------------------------------------
modulus.py
# routines for modulo calculations
# written by Jurjen Bos (jurjen@q2c.nl)
# uses random.py for generation of random numbers
# defines class modular that describes a modulus for computations
# and class modularNumber that contains a number with modulus
# the chinese remainder theorem is used. This implies that if the prime
# factorization of the modulus is known, functions will be fast, and that
# several extra functions (e.g., computation of roots) are allowed.
# the modulus is stored as a list of moduli that are multiplied to get the
modulus
# functions fracPower, phi of modular and root of modularNumber assume
# that all moduli are prime
import modulus
import random
def product(s): #compute product of elements of sequence
result = 1
for i in s: result = result*i
return result
def divMod(a, b, n): # division modulo n; From Knuth, Volume 2,
page 325 (u is b, v is n)
u1, u3, v1, v3 = a, b, 0, n
while v3<>0:
q = u3/v3
(u1, u3), (v1, v3) = (v1, v3), (u1-v1*q, u3-v3*q)
if u3 <> 1: raise ZeroDivisionError, 'modular division'
return u1%n
def powMod(a, b, n): # power modulo n
result = 1
while b>0:
while b%2==0:
a, b = a*a%n, b/2
result, b = result*a%n, b-1
return result
def computeUnits(n): # compute values of unit vectors in tuple
representation
p = product(n)
result = []
for d in n: result.append(p/d*divMod(1, p/d%d, d))
return result
class modular:
def __init__(self, *n):
self.moduli = n
self.product = product(n)
self.units = computeUnits(n)
self.rootCache = {}
def __repr__(self):
if len(self.moduli)==1: return 'modular(' +
`self.moduli[0]` + ')'
return 'modular' + `self.moduli`
def number(self, a): # create modularNumber with value a
return modularNumber(self, a)
def random(self): return modularNumber(self, random.rand(self.product))
def split(self, a): # split value in moduli
s = []
for d in self.moduli: s.append(a%d)
return s
def combine(self, s): # combine moduli to value (chinese
remainder theorem)
a = 0
if len(self.units)<>len(s): raise ValueError, 'number of
elements'
for i in range(len(self.units)): a = a+self.units[i]*s[i]
return a%self.product
def fracPower(self, n, d): # convert n/d power to powers for
each modulus
# assumes moduli are prime!
if self.rootCache.has_key((n,d)): return self.rootCache[(n,d)]
result = []
for m in self.moduli: result.append(divMod(n, d, m-1))
self.rootCache[(n,d)] = result
return result
def phi(self): # euler's phi function.
Assumes moduli are prime!
result = 1
for m in self.moduli: result = result * (m-1)
return result
class modularNumber(modular):
def __init__(self, modulus, s):
self.moduli = modulus.moduli
self.product = modulus.product
self.units = modulus.units
self.rootCache = modulus.rootCache
if type(s) == listType: self.values = s
else: self.values = self.split(s)
if len(self.values)<>len(self.moduli): raise ValueError,
'number of elements'
def __repr__(self):
return `self.value()` + ' (mod ' + `self.product` + ')'
def power(self, d): # power. if moduli are not prime,
only works if 0<=d<modulus
s = []
for i in range(len(self.moduli)):
s.append(powMod(self.values[i],
d%(self.moduli[i]-1), self.moduli[i]))
return modularNumber(self, s)
def root(self, d): # root
s = []
r = self.fracPower(1, d)
for i in range(len(self.moduli)):
s.append(powMod(self.values[i], r[i], self.moduli[i]))
return modularNumber(self, s)
def value(self): # return the value
return self.combine(self.values)
def __add__(self, other):
if type(other)<>modType:
return self+modularNumber(self, other)
if other.moduli is not self.moduli: raise ValueError, '+
with differing modulus'
s = []
for i in range(len(self.moduli)):
s.append((self.values[i]+other.values[i])%self.moduli[i])
return modularNumber(self, s)
def __sub__(self, other):
if other.moduli is not self.moduli: raise ValueError, '-
with differing modulus'
s = []
for i in range(len(self.moduli)):
s.append((self.values[i]-other.values[i])%self.moduli[i])
return modularNumber(self, s)
def __mul__(self, other):
if type(other)<>modType:
return self*modularNumber(self, other)
if other.moduli is not self.moduli: raise ValueError, '*
with differing modulus'
s = []
for i in range(len(self.moduli)):
s.append((self.values[i]*other.values[i])%self.moduli[i])
return modularNumber(self, s)
def __div__(self, other):
if other.moduli is not self.moduli: raise ValueError, '/
with differing modulus'
s = []
for i in range(len(self.moduli)):
s.append(divMod(self.values[i], other.values[i],
self.moduli[i]))
return modularNumber(self, s)
def __coerce__(self, other):
return self, modularNumber(self, other)
def __pos__(self): return self
def __neg__(self):
s = []
for i in range(len(self.moduli)):
s.append((-self.values[i])%self.moduli[i])
return modularNumber(self, s)
def __cmp__(self, other): # note that < and > have no really
useful meaning
for i in range(len(self.moduli)):
if self.values[i]<>other.values[i]: return
cmp(self.values[i], other.values[i])
return 0
def __hash__(self):
return hash(self.values)^hash(self.moduli)
def __nonzero__(self):
for a in self.values:
if a: return 1
return 0
def __hex__(self): return hex(self.value())
listType = type([])
modType = type(modularNumber(modular(), 0))
/*--- Jurjen N.E. Bos ---------------------- Email: jurjen@q2c.nl ---*/
int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;for(;d=0,g=c*2;c
-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}