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