Calculate Pi with Python

Intro to Computer programming worked at calculating digits of pi today. The actual algorithms aren’t too bad, but getting more than the standard number of digits from a double is a bit trickier. Here’s a program that calculates pi using:

Bailey–Borwein–Plouffe formula

Bellard’s formula

and
Chudnovsky algorithm

Holy smokes is Chudnovsky algorithm’s fast!

			 Plouff 		 Bellard 			 Chudnovsky
Iteration number  1   3.133333333333333333333333   3.141765873015873015873017   3.141592653589734207668453
Iteration number  2   3.141422466422466422466422   3.141592571868390306374053   3.141592653589793238462642
Iteration number  3   3.141587390346581523052111   3.141592653642050769944284   3.141592653589793238462642
Iteration number  4   3.141592457567435381837004   3.141592653589755368080514   3.141592653589793238462642
Iteration number  5   3.141592645460336319557021   3.141592653589793267843377   3.141592653589793238462642
Iteration number  6   3.141592653228087534734378   3.141592653589793238438852   3.141592653589793238462642
Iteration number  7   3.141592653572880827785241   3.141592653589793238462664   3.141592653589793238462642
Iteration number  8   3.141592653588972704940778   3.141592653589793238462644   3.141592653589793238462642
Iteration number  9   3.141592653589752275236178   3.141592653589793238462644   3.141592653589793238462642
Iteration number  10   3.141592653589791146388777   3.141592653589793238462644   3.141592653589793238462642
Iteration number  11   3.141592653589793129614171   3.141592653589793238462644   3.141592653589793238462642
Iteration number  12   3.141592653589793232711293   3.141592653589793238462644   3.141592653589793238462642
Iteration number  13   3.141592653589793238154767   3.141592653589793238462644   3.141592653589793238462642
Iteration number  14   3.141592653589793238445978   3.141592653589793238462644   3.141592653589793238462642
Iteration number  15   3.141592653589793238461733   3.141592653589793238462644   3.141592653589793238462642
Iteration number  16   3.141592653589793238462594   3.141592653589793238462644   3.141592653589793238462642
Iteration number  17   3.141592653589793238462641   3.141592653589793238462644   3.141592653589793238462642
Iteration number  18   3.141592653589793238462644   3.141592653589793238462644   3.141592653589793238462642
Iteration number  19   3.141592653589793238462644   3.141592653589793238462644   3.141592653589793238462642

Source code (Python)

from decimal import *

#Sets decimal to 25 digits of precision
getcontext().prec = 25

def factorial(n):
    if n<1:
        return 1
    else:
        return n * factorial(n-1)

def plouffBig(n): #http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(1)/(16**k))*((Decimal(4)/(8*k+1))-(Decimal(2)/(8*k+4))-(Decimal(1)/(8*k+5))-(Decimal(1)/(8*k+6)))
        k += 1
    return pi

def bellardBig(n): #http://en.wikipedia.org/wiki/Bellard%27s_formula
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1) + Decimal(1)/(10*k+9) - Decimal(64)/(10*k+3) - Decimal(32)/(4*k+1) - Decimal(4)/(10*k+5) - Decimal(4)/(10*k+7) -Decimal(1)/(4*k+3))
        k += 1
    pi = pi * 1/(2**6)
    return pi

def chudnovskyBig(n): #http://en.wikipedia.org/wiki/Chudnovsky_algorithm
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k)))
        k += 1
    pi = pi * Decimal(10005).sqrt()/4270934400
    pi = pi**(-1)
    return pi
print "\t\t\t Plouff \t\t Bellard \t\t\t Chudnovsky"
for i in xrange(1,20):
    print "Iteration number ",i, " ", plouffBig(i), " " , bellardBig(i)," ", chudnovskyBig(i)


Posted in Full Posts | Tagged , , , | 9 Comments

Powers of 2 without a Zero


Interesting. The Intro to Programming Students gave a crack at this one; and after 30 minutes all were close and half had it. Some worked in Python (easier large numbers and string manipulation) and some in Java (muuuuch faster). Here’s my Java (BlueJ) code. It checks all powers of 2 up to 2^10,000, and then it prints out the counts of all the digits 0-9. The digits are spread out pretty evenly.

import java.math.BigInteger;
public class Powersof2
{
    public int[] digits = new int[11];;
    public Powersof2()
    {
        BigInteger number = new BigInteger ("1");;
        int power = 1;
        
        for (int i = 0; i < digits.length; i++)
        {
            digits[i]=0;
        }
        
        while (power <= 10000)
        {
            number = number.multiply(new BigInteger ("2"));
            
            if (HasAZeroInIt(number))
            {
                System.out.println("Next power of 2 is " + number + " with an exponent of " + power);
            }
            
            
            power += 1;
        }
        
        for (int i = 0; i < digits.length-1; i++)
        {
            System.out.println("index " + i + " and number of digits is " + digits[i]);
        }
                   
       
    }
    
    public boolean HasAZeroInIt(BigInteger n)
    {
        String strn = String.valueOf(n);
        
        int index = 0;
        while (index < strn.length())
        {
            digits[Character.getNumericValue(strn.charAt(index))] += 1;
            index += 1;
        }
        index = 0;
        while (index < strn.length())
        {
            if (strn.charAt(index) == '0')
            {
                return false;
            }
            index += 1;
        }
        return true;
    }    
}
Posted in Uncategorized | Tagged , , , | Leave a comment

Reducing Prime

Otherwise known as left-truncating primes with every suffix prime and no digit are zero. Here is the whole sequence. It was proven that 357686312646216567629137 is the last number of the sequence in 1999. Cool that math understandable to me is still being done in my lifetime. Although I’d bet that the actual proof is beyond my pay grade.

Here’s a recursive python program to print out a list of these numbers under 1 million. Unfortunately it doesn’t work on the number in question, my prime checker isn’t efficient enough.

from math import sqrt
#function to check to see if a number is prime
def prime(n):
    for i in xrange(2,int(sqrt(n))+1):
        if (n%i == 0):
            return False
        return True

#recursive function to check to see if a number continues to be prime
#as the left digit gets taken away
def reducingprime(n):
    #Stopping condition
    if n == 0:
        return True
    else:
        string_n = str(n)
        if n < 10:
            return prime(n)

        #check recursive step only if second digit isn't 0
        elif string_n[1]!='0' and prime(n):
            #this takes away the left digit
            next = n % (10**((len(str(n))-1)))
            return reducingprime(next)
        else:
            return False


for a in xrange(2,1000000):
    if reducingprime(a):
        while a > 0:
            print a," ",
            a = a%(10**((len(str(a))-1)))
        print " "
Posted in interesting stuff | Tagged , , | Leave a comment

Oreo: Original vs. Double vs. Mega

[edit 8/20/13: New post with more data on the weights of Oreos and Double Stuf Oreos.]
[edit 8/21/13: And Part 3.]

Are Double Stuf Oreos really double stuffed? Some disagreement here: Chris Danielson in part 1 2 and 3 versus Chris Lusto in part 1 and 2.

Enter in the Mega Stuf Oreo.
megastuf

Is the Mega Stuf double stuffed, triple stuffed, more???? And who was correct, Chris or Chris?

Evidence:

10 Original Oreos:
Original Oreos

10 Double Stuf Oreos:
Double Stuf

10 Mega Stuf Oreos:
Mega Stuf

5 Wafers:
Wafers

I’ll leave the math to you. Or if you’re terrifically lazy, spoilers found here.

Edit (8/18/2013): This post is now featured on the Huffington Post. Should have used better handwriting!

Posted in Full Posts | Tagged , , | 93 Comments

A Billion Heartbeats

Interesting.

As animals get bigger, from tiny shrew to huge blue whale, pulse rates slow down and life spans stretch out longer, conspiring so that the number of heartbeats during an average stay on Earth tends to be roughly the same, around a billion.

And …

… the surface area a creature uses to dissipate the heat of the metabolic fires does not grow as fast as its body mass.

To see this, consider a mouse as an approximation of a small sphere. As the sphere grows larger, to cat size, the surface area increases along two dimensions but the volume increases along three dimensions. The size of the biological radiator cannot possibly keep up with the size of the metabolic engine.

Haha. Ass.
Might be fun to start with the theory that all animals get a billion heartbeats and build up some charts to see if its true. Also a surface area of various mammals vs age graph might be interesting to calculate and plot.
Does this mean that I should stop exercising?
Link from Kottke.

Posted in interesting stuff | Tagged | 2 Comments

Spirograph

Fantastic.

Questions for students:

  • What’s with that fraction on the left side? What does that do?
  • How do I get a graph with 7 “arms”?
  • What is a the graph with 2 arms called?
  • How many circuits do I need to finish the shape if my fraction is 25/60?

 

Posted in interesting stuff | Tagged , | 1 Comment