I was inspired by the Pi program that is distibuted with Python. I decided to expand it a
little, and speed it up using a better approximation method.
The result is a program that can print the following number to infinite precision:
- pi
- e to fraction powers
- sin, cos, sinh, cosh, asin, atan, acos of a fraction
- roots of fractions
- sums of any of the above
At the end of the program, there are some examples. For example, if you are only
interested in printing pi, use p(pi4).
Hope you like it; let me know if you do.
# This program prints mathematical constants to infinite precision
# It can print all constants specifiable as num/den*sum(i=0,Β,product(k=1,i,f(n)/g(n)))
# where f and g integral functions satisfying -g(k),
# or sums of these constants.
# Written by Jurjen N.E. Bos (jurjen(a)q2c.nl), 1996-97.
# Examples are shown below
from sys import stdout
from string import zfill
# the class below is dirty in the sense that class instances modify themselves
# after quite "innocent" operations, such as *,/,+
negDenom = "Denominator must be positive"
class series:
def __init__(self,num,den,f,g):
self.num = long(num) # current value (always positive)
self.term = self.num # next term (may be negative)
self.den = long(den) # common denominator (always positive)
self.f = f
self.g = g
self.k = 1
# multiply sequence by a number (modifying the sequence!)
def __mul__(self,n):
self.num = n*self.num
self.term = n*self.term
return self
# divide sequence by a number (modifying the sequence!)
def __div__(self,n):
self.den = n*self.den
return self
# add number to sequence (modifying the sequence!)
def __add__(self,n):
self.num = self.num+n*self.den
return self
# expand with given precision
def expand(self,factor):
while abs(self.term)*factor>=self.den:
self.term = self.term*self.f(self.k)
g=self.g(self.k)
self.num = self.num*g
self.den = self.den*g
self.num = self.num+self.term
self.k = self.k+1
# this is a lower estimate: subtract one if negative term
return self.num*factor/self.den-(self.term<0)
# integer part of sum of a list of series
def intPart(l):
extra = 1
# force the while to execute once
estimate = 0
# we have fractional part and error for every term,
# so we can expect up to 2 extra per term
while estimate%extra >= extra-3*len(l):
# too much error, estimate two more digits
extra = extra*100L
estimate = 0
for s in l:
estimate = estimate+s.expand(extra)
return int(estimate/extra)
def p(l):
digits = 8 # number of digits to print at once
base = 10**digits # our computation "number base"
lastInt = intPart(l)
stdout.write(`lastInt`+',')
#print fraction
while 1:
carry = lastInt
for s in l:
q,s.num = divmod(s.num,s.den)
carry = carry-int(q)
s*base
lastInt = intPart(l)
output = lastInt-carry*base
if not (0<=output<base): raise "Consistency check failed"
stdout.write(zfill(output,digits))
stdout.flush()
# now some interesting series
# e to the power n/m
def expSeries(n=1,m=1):
if m<0: raise negDenom
def f(k,nn=n): return nn
def g(k,mm=m): return k*mm
return series(1,1,f,g)
# c*arc sine of n/m
def asinSeries(n,m=1):
if m<0: raise negDenom
def f(k,n2=n**2): return (2*k-1)**2*n2
def g(k,m2=m**2): return m2*(2*k)*(2*k+1)
return series(n,m,f,g)
# arc tangent of 1/n
def atanSeries(n,m=1):
if m<0: raise negDenom
def f(k,n2=n**2): return -(2*k-1)*n2
def g(k,m2=m**2): return m2*(2*k+1)
return series(n,m,f,g)
# r-th root of n/m (around one)
def rootSeries(n,m,r=2):
if m<0: raise negDenom
n = n-m # move origin to zero
def f(k,nn=n,rr=r): return nn*(1-(k-1)*rr)
def g(k,mr=m*r): return k*mr
return series(1,1,f,g)
# the sine of n/m
def sinSeries(n,m=1):
if m<0: raise negDenom
def f(k,n2=n**2): return -n2
def g(k,m2=m**2): return (2*k)*(2*k+1)*m2
return series(n,m,f,g)
# the cosine of n/m
def cosSeries(n,m=1):
def f(k,n2=n**2): return -n2
def g(k,m2=m**2): return (2*k-1)*(2*k)*m2
return series(1,1,f,g)
# the hyperbolic sine of n/m
def sinhSeries(n,m=1):
if m<0: raise negDenom
def f(k,n2=n**2): return n2
def g(k,m2=m**2): return (2*k)*(2*k+1)*m2
return series(n,m,f,g)
# the hyperbolic cosine of n/m
def coshSeries(n,m=1):
def f(k,n2=n**2): return n2
def g(k,m2=m**2): return (2*k-1)*(2*k)*m2
return series(1,1,f,g)
# natural log of n/m (around one)
def logSeries(n,m):
if m<0: raise negDenom
n = n-m # move origin to zero
def f(k,nn=n): return -k*nn
def g(k,mm=m): return mm*(k+1)
return series(n,m,f,g)
# now try a few of these: just type p(<value>), e.g. p(pi4)
# (this works only once for each <value>, since it is modified!)
# the constant e is the simplest
e=[expSeries()]
# pi using some standard formulas (in speed order: last one is fastest)
pi1=[series(4,2,lambda(k): k,lambda(k): 2*k+1)]
pi2=[asinSeries(1,2)*6]
pi3=[atanSeries(1,5)*16,atanSeries(1,239)*-4]
pi4=[atanSeries(1,8)*24,atanSeries(1,57)*8,atanSeries(1,239)*4]
# roots are easily speeded up by rational approximations of the result
# square root of 2
sr2=[rootSeries(2*14**2,17**2)*17/12]
# cube root of 2
cr2=[rootSeries(2*4**3,5**3,3)*5/4]
# golden ratio
phi=[(rootSeries(5*17**2,38**2)*38/17+1)/2]
# log of 2 (simple way)
ln2a=[logSeries(1,2)*-1]
# log of 2,3,5,7 (fast way, using combinations of small logs)
def combineLogs(a,b,c,d):
return [logSeries(125,126)*-a,logSeries(224,225)*-b,
logSeries(2400,2401)*c,logSeries(4374,4375)*-d]
ln2=combineLogs( 72,27,19,31)
ln3=combineLogs(114,43,30,49)
ln5=combineLogs(167,63,44,72)
ln7=combineLogs(202,76,53,87)
# log of 10 is sum of log 2 and 5
ln10=combineLogs(239,90,63,103)
-- jurjen(a)q2c.nl ( Jurjen N.E. Bos)
========PYTHONMAC-SIG - SIG on Python for the Apple Macintosh
send messages to: pythonmac-sig(a)python.org
administrivia to: pythonmac-sig-request(a)python.org
========