Power Pi


I was asked recently how people calculate PI. That is, what method is employed on a supercomputer to get the value of PI to a specific number of decimal places. I didn’t know off-hand, so I looked it up. Wolfram gives 32 power expansions for PI. For those of you unfamiliar with the math, a power expansion (or power series) is just another way to write out a number as the summation of individual arithmetic operations (i.e. – terms).

The simplest one to implement in VBScript is:
pi = 4 * sum [(-1)^k / (2*k + 1)] where k goes from 0 to infinity in steps of 1.

The (-1)^k part causes each iteration of the formula to oscillate between two points that gradually close in on the correct value. Then the 1/(2*k+1) is the portion of the power expansion that actually contributes to the final value. It’s easy to write, but it is slow. On my laptop, I need to run k up to 5,000,000 terms, and it takes 15 seconds to finish. Even then, the result is only accurate to 6 decimal places.

I tried to streamline the program a bit, so the part where I multiply the result by 4 isn’t until the final wscript.echo statement. Note also that there’s little point in calculating (-1)^k EVERY SINGLE TIME, so I just use sn2 to hold the current value of the sign flipping component.

Lastly,
for k=0 to 500
tmp = sn2/(2*k + 1)
next

is the same as

for k=1 to 1001 Step 2
tmp = sn2/k
next

but is maybe 25% faster to run.

res = 0
sn = -1
sn2 = 1
max = 2 * 5000000
for k = 1 to max Step 2
tmp = sn2 / k
res = res + tmp
sn2 = sn2 * sn
next

wscript.echo “res = ” & cstr(4 * res) & ” term ” & cstr(k) & ” = +/-” & cstr(2 * tmp)

The actual value of PI is: 3.14159265359
My value is 3.14159285359 At step 5,000,000
and 3.14159245359 At step 5,000,001
So, I’m at +/- 0.000000199999.

Not a really useful program, but it is fun to see the answer coming out close.

—-

There’s a more direct method with just two terms (using Excel expressions):

pi = 4 * (4 * atan(1/5) – atan(1/239))

This of course now involves finding the irrational values of two separate numbers to whatever accuracy you’re seeking for pi… Maybe not an improvement.

—-

Many of the formulas described at the Mathworld page are circular, in that to find pi, you have to use pi within the power series. Which to me seems to be kind of hard to do in software. However, I did like the version which comes from the special value of the Riemann zeta function for (pi^2)/6 = zeta(2n) for all positive integers n:

pi = sqrt(6*(1 + 1/2^2 + 1/3^2 + 1/4^2 + …))

In VBScript, this becomes:

t = 0
for n = 1 to 5000000
t = t + 1/n^2
next
wscript.echo “pi = ” & cstr((6 * t)^0.5)

pi = 3.14159246260382,

This is 29% faster than the first program above (5 seconds for 500 million terms, compared to 6.5 seconds for the first program), and they are both still accurate to only 6 places. (If both programs run for 500 million terms, the first program takes 357.6 seconds, and the second for 255.7 seconds, giving program #2 a 28.6% speed advantage, although #2 was only accurate to 7 places and #1 was already accurate to 8 places.)

After 500 million terms:
PI = 3.14159265359
#1 = 3.14159265158926
#2 = 3.14159264498239

Advertisements
Leave a comment

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: