# Three useful classes

Jurjen N.E. Bos (jurjen@q2c.nl)
Mon, 28 Mar 1994 11:13:08 +0200

Hi,
Guido told me that the python community might be happy with my homebrew
classes, so I publish them here for your enjoyment. You may use them
freely under the following conditions:
- If you have improvements or comments, let me know!
- Leave my name in it: maybe it will make me rich and famous :-)

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

#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
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
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)

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);}