My first Python program.. Calculating pi geometrically.

[Ron Adam]

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:  radam2 at tampabay.rr.com or ron at 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 #