code optimization (calc PI)
Nick Craig-Wood
nick at craig-wood.com
Thu Jan 4 07:30:07 EST 2007
mm <michael at mustun.ch> wrote:
> (Yes, I konw whats an object is...)
> BTW. I did a translation of a pi callculation programm in C to Python.
> (Do it by your own... ;-)
> Calc PI for 800 digs(?). (german: Stellen)
> 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);}
Below you'll find my translation. Note that I improved the algorithm
greatly. Note the nice easy to read algorithm also! This calculated
800 digits of pi in 0.1 seconds and 10,000 digits in 20 seconds.
------------------------------------------------------------
"""
Standalone Program to calculate PI using python only
Nick Craig-Wood <nick at craig-wood.com>
"""
import sys
from time import time
class FixedPoint(object):
"""A minimal immutable fixed point number class"""
__slots__ = ['value'] # __weakref__
digits = 25 # default number of digits
base = 10**digits
def __init__(self, value=0):
"""Initialise the FixedPoint object from a numeric type"""
try:
self.value = long(value * self.base)
except (TypeError, ValueError), e:
raise TypeError("Can't convert %r into long", value)
def set_digits(cls, digits):
"""Set the default number of digits for new FixedPoint numbers"""
cls.digits = digits
cls.base = 10**digits
set_digits = classmethod(set_digits)
def copy(self):
"Copy self into a new object"
new = FixedPoint.__new__(FixedPoint)
new.value = self.value
return new
__copy__ = __deepcopy__ = copy
def __add__(self, other):
"""Add self to other returning the result in a new object"""
new = self.copy()
new.value = self.value + other.value
return new
__radd__ = __add__
def __sub__(self, other):
"""Subtract self from other returning the result in a new object"""
new = self.copy()
new.value = self.value - other.value
return new
def __rsub__(self, other):
new = self.copy()
new.value = other.value - self.value
return new
def __mul__(self, other):
"""
Multiply self by other returning the result in a new object. The
number of digits is that of self.
Note that FixedPoint(x) * int(y) is much more efficient than
FixedPoint(x) * FixedPoint(y)
"""
new = self.copy()
if isinstance(other, (int, long)):
# Multiply by scalar
new.value = self.value * other
else:
new.value = (self.value * other.value) // other.base
return new
__rmul__ = __mul__
def __div__(self, other):
"""
Divide self by other returning the result in a new object. The
number of digits is that of self.
Note that FixedPoint(x) / int(y) is much more efficient than
FixedPoint(x) / FixedPoint(y)
"""
new = self.copy()
if isinstance(other, (int, long)):
# Divide by scalar
new.value = self.value // other
else:
new.value = (self.value * other.base) // other.value
return new
def __rdiv__(self, other):
if not isinstance(other, FixedPoint):
other = FixedPoint(other, self.digits)
return other.__div__(self)
def __str__(self):
"""
Convert self into a string. This will be the integral part of
the number possibly preceded by a - sign followed by a . then
exactly n digits where n is the number of digits specified on
creation.
"""
if self.value < 0:
s = str(-self.value)
else:
s = str(self.value)
s = s.zfill(self.digits + 1)
s = s[:-self.digits] + "." + s[-self.digits:]
if self.value < 0:
s = "-" + s
return s
def __abs__(self):
"""Return the absolute value of self"""
new = self.copy()
if new.value < 0:
new.value = - new.value
return new
def __cmp__(self, other):
"""Compare self with other"""
return cmp(self.value, other.value)
def sqrt(n):
"""
Return the square root of n. It uses a second order
Newton-Raphson convgence. This doubles the number of significant
figures on each iteration. This will work for any finite
precision numeric type.
"""
x = n / 2 # FIXME find a better guess!
old_delta = None
while 1:
x_old = x
x = (x + n / x) / 2
delta = abs(x - x_old)
if old_delta is not None and delta >= old_delta:
break
old_delta = delta
return x
def pi_agm(one = 1.0, sqrt=sqrt):
"""
Brent-Salamin Arithmetic-Geometric Mean formula for pi. This
converges quadratically.
FIXME can increase the number of digits in each iteration to speed up?
"""
_1 = one
_2 = _1 + _1
a = _1
b = _1/sqrt(_2)
t = _1/2
k = 1
two_to_the_k = 2L
old_delta = None
while 1:
s = (a + b ) / 2
d = a - s
d = d * d
a = s
s = s * s
t = t - two_to_the_k * d
b = sqrt(s - d)
delta = abs(a - b)
# print k, delta, old_delta
if old_delta is not None and delta >= old_delta:
break
old_delta = delta
k = k + 1
two_to_the_k *= 2
a = a + b
return ( a * a ) / ( _2 * t)
if __name__ == "__main__":
try:
digits = int(sys.argv[1])
except:
digits = 100
print "Calculating",digits,"digits of pi"
start = time()
FixedPoint.set_digits(digits)
_1 = FixedPoint(1)
print pi_agm(_1, sqrt=sqrt)
dt = time() - start
print "That took",dt,"seconds"
------------------------------------------------------------
> But Python is much slower here then C here. I used a while-loop within a
> while-loop. Not that much calculation here.
There is more than one way to skin this cat. A better algorithm helps
a lot.
If I really wanted lots of digits of pi from python I'd do it like
this. This is an example of the right tool for the job. 0.5ms looks
pretty good to me!
>>> import math, gmpy
>>> bits = int(800/math.log10(2))
>>> start = time(); pi=gmpy.pi(bits); dt = time()-start
>>> print "That took",dt,"seconds"
That took 0.00055193901062 seconds
>>> pi
mpf('3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455e0',2657)
>>>
--
Nick Craig-Wood <nick at craig-wood.com> -- http://www.craig-wood.com/nick
More information about the Python-list
mailing list