Hi, I'm having fun learning Python and want to say thanks to everyone
here for a great programming language.
Below is my first Python program (not my first program) and I'd
apreciate any feedback on how I might do things differently to make it
either more consice, readable, or faster. ie... are there better ways
to do it in Python?
It won't break any records for calculating pi, that wasn't my goal,
learning Python was. But it might be good for testing floating point
and decimal performance. It also might make a good demo of the same.
Anyway, here it is, any feedback would be apreciated.
_Ron
"""
## Program name: pi10.py
## Author: Ronald R. Adam
## Version: 1.01
## Date: 11/01/2003
##
## Email:
ra****@tampabay.rr.com or
ro*@ronadam.com
##
## This program calculates Pi using a geometric progression of
## triangles. The program is written in Python and requires
## the gmpy math module in order to use extended floating point
## precision.
##
## This is a rough diagram of the location of the points used.
##
## |
## A* --.__
## | \ `-._
## | \ `* D
## | \ / \
## | C * \
## | . \ \
## | . \ |
## | . \ |
## *---------------*---
## O B
##
## Pint C is the mid point of line AB.
## Point D is on the circle 90 degrees to line AB from point C.
##
## Area = OAB*4 + CDB*8 + ...
##
## ... Point A is relocated to point D and then new points
## C and D are calculated. This is repeated ...
##
## ... + CDB*16 + CDB*32 + ... (Repeats until done)
##
## When it's done, the last line of AB length multiplied by with
## the number of sides (2^n) is used to calculate the
## circumference, and then calculate pi again to check the result.
##
## Pi to 10,000 digits take about 5 minutes on my desktop computer.
## I'm not sure what the maximum is so drop me an email if you get
## some really long sequences. This program progrssivly slows down
## as the number of digits increases because it uses the entire
## length of each floating point number.
"""
import gmpy
import datetime
import time
print """
Pi10.py Version 1.0, by Ronald R. Adam
This program calculates Pi by using a geometric progression
of triangles to measure the area of a circle with a radius of 1.
"""
d=1
while d>0:
try:
d = gmpy.mpz(raw_input('Number of decimal digits: '))
except:
break
# Estimate the precision needed.
# This may need some adjustment.
p = int(d*3.325)+16 # gmpy precision value
start = time.clock()
t_old = start
print 'Calculating Pi ...'
# set variable values and precision.
n = gmpy.mpz('2') # number of triangles
radius = gmpy.mpf('1.0',p) # radius
area = gmpy.mpf((radius*radius)/2 * 2**n) #area=(base*height)/2*4
Ax, Ay = gmpy.mpf('0.0',p) , radius # point A
Bx, By = radius , gmpy.mpz('0') # point B
Cx = Cy = gmpy.mpf('0.0',p) # point C
OC = CD = CB = gmpy.mpf('0.0',p) # lines OC,CD,CB
vx = vy = gmpy.mpf('0.0',p) # vectors of line CD
# Loop until additions to area are too small to calculate
# using the requested precision + 2 digits.
p_area = gmpy.mpf('0.0',p)
while gmpy.fdigits(p_area,10,d+2) <> gmpy.fdigits(area,10,d+2):
p_area = area
Cx, Cy = (Ax+1)/2, Ay/2 # find midpoint of line AB
OC = gmpy.fsqrt(Cx**2 + Cy**2) # find length of line OC
CD = radius - OC # find length of line CD
CB = gmpy.fsqrt((Bx-Cx)**2+Cy**2) # find length of line CB
n += 1 # double triangles 2^(n+1)
area = area + (CB*CD/2) * 2**n # add trianges BCD * 2^n
vx = CD * Cx/OC # find x vector of line CD
vy = CD * Cy/OC # find y vector of line CD
Ax, Ay = Cx+vx, Cy+vy # Ax,Ay = Dx,Dy
t = time.clock()
if t > t_old+10:
t_old = t
print 'Exceeding 2^', n ,'sides ...'
# pi1 = area/radius^2
pi1 = gmpy.fdigits(area) # use (area/radius**2) if radius <>1
# pi2 = circumference /(2*radius)
pi2=gmpy.fdigits((gmpy.fsqrt((Bx-Ax)**2+(By-Ay)**2)*2**n)/(2*radius))
# Compare the two above calculations to each other
# and copy the matching & requested digits to the
# result string.
pi = ''
k = 0
while pi1[k] == pi2[k]:
pi = pi+pi1[k]
k += 1
pi = pi[0:int(d+2)]
# Report the result
print
print 'Pi =\n',pi
print
print len(pi)-2, 'decimal places'
print 'Radius =', radius
print '2^', n, ' sides'
stop = time.clock()
print 'Elapsed time:',datetime.timedelta( 0, stop-start)
print
# End #
