nth digit of PI

BlitzMax Forums/BlitzMax Beginners Area/nth digit of PI

Shortwind(Posted 2010) [#1]
Hi, I'm trying to program the formula to find the Nth digit of PI.

Here is one version of the formula, but I'm doing something wrong and can't make it work...

Quoted from http://www.math.hmc.edu/funfacts/ffiles/20010.5.shtml
Here is a very interesting formula for pi, discovered by David Bailey, Peter Borwein, and Simon Plouffe in 1995:

Pi = SUMk=0 to infinity 16-k [ 4/(8k+1) - 2/(8k+4) - 1/(8k+5) - 1/(8k+6) ].

...

Now Wolfram says something different:

http://mathworld.wolfram.com/PiDigits.html

Xn = frac( 16Xn-1 + 120n^2 - 89n +16 / 512n^4 - 1024n^3 + 712n^2 - 206n - 21)

I can't get either of these to work. It's got to be my interpretation of the formulas.

My non-working attempt...


'4/(8k+1) - 2/(8k+4) - 1/(8k+5) - 1/(8k+6)
Local k:Double
Local j:Double
Local x:Double

For k=0 To 4

j=16^-k*( (4/(8*k+1)) - (2/(8*k+4)) - (1/(8*k+5)) - (1/(8*k+6)))

Print j

Next 

'Xn = frac( 16Xn-1 + (120n^2 - 89n +16 / 512n^4 - 1024n^3 + 712n^2 - 206n - 21))
For k=0 To 4

x= 16*(k-1) + (((120*n^2) - (89*n +16)) / ((512*n^4) - (1024*n^3) + (712*n^2) - (206*n) - 21))

Print x

Next



skidracer(Posted 2010) [#2]
you forgot to sum the first one:

'4/(8k+1) - 2/(8k+4) - 1/(8k+5) - 1/(8k+6)
Local k:Double
Local j:Double
Local x:Double

For k=0 To 24

j:+16^-k*( (4/(8*k+1)) - (2/(8*k+4)) - (1/(8*k+5)) - (1/(8*k+6)))

Print j

Next 



Shortwind(Posted 2010) [#3]
Thanks skidracer, but I guess I'm still confused.

If I want to find the 29th digit of PI, I should be able to make k=29, send it to the formula and have 7 as the answer. Do I need to change to a different print routine to make it work correctly?

The whole point of this formula is to make it so you don't have to calculate all the previous digits. I just stuck it in a loop to see what answers I'd get.

Maybe I'm just not smart enough to figure this out.... ho hummmm.

:(


Floyd(Posted 2010) [#4]
There's a lot of detail missing from that description. And it ultimately gives base 16 digits.

Pi in base 10 is 3.14159..., which means 3 + 1/10 + 4/100 + 1/1000 etc.

Pi in base 16 is 3.243F6..., which is 3 + 2/16 + 4/256 + 3/4096 + 15/65536 etc.

First, here's a demonstration that the infinite sum gives the value of Pi.
Then we pick out the first few hexadecimal digits.

Local k:Double
Local j:Double
Local x:Double
Local Pie:Double

For k=0 To 15  ' not quite infinity, but close enough. 
	Pie :+ 16^-k*( (4/(8*k+1)) - (2/(8*k+4)) - (1/(8*k+5)) - (1/(8*k+6)))
Next 

Print
Print " 0 to 15 sum: " + Pie
Print " Value of Pi: " + Pi
Print 

' Now let's see some hexadecimal digits of Pi.
Local n:Int, d:Int, t:String

For n = 0 To 10
	d = Int(Pie)
	Pie :- d
	Pie :* 16
	t :+ "0123456789ABCDEF"[d..d+1]
Next
Print "Pi in base 16: " + t[0..1] + "." + t[1..] + "..."


The point of all this is that the initial sum to compute Pi contains powers of 1/16. At first it looks like you need to do the entire sum. That's what I do in the sample code above. But with careful estimates you can figure out what a particular digit is with much less computation. As I said, there are a lot of details to fill in.

Also, you really can't get very far with double precision floats. They give 53-bit accuracy. Calculating Pi out to the 29th hexadecimal digit needs far more as each digit is four bits. Allowing for some extra precision to ensure rounding to the proper digit you would probably need about 130-bit accuracy.

Last edited 2010


Floyd(Posted 2010) [#5]
This is interesting. It looks like that MathWorld formula has a typo.

I implemented it and the resulting digits looked nothing like reality.

Notice that the terms in the numerator go + - +. The denominator goes + - + - -.
We can reasonably speculate that should be + - + - +.

Now the result agrees with the code I posted earlier: 3.243F6A8885
Local x:Double, n:Double, d:Int

For n = 1 To 10	
	x = 16*x + (120*n^2 - 89*n +16) / (512*n^4 - 1024*n^3 + 712*n^2 - 206*n + 21)
	x :- Floor(x)   ' x = Frac(x)
	d = Int(16*x)   ' nth hexadecimal digit after 3.
	Print d
Next
Note this can't be continued for long. Every time through the loop we multiply by 16 and discard the integer part. This throws away 4 bits. So with 53-bit arithmetic we can get twelve digits, maybe thirteen. After that the results are certainly junk.