Methods of Estimating Pi

[Orion1]

I agree, the fault is mine for selecting the Maclaurin series for arctan(x) for my initial evaluation, and that further discussion on this series needs to be suppressed because it represents a massive sidetrack.

 

 

[math]\pi = \sum_{n=0}^{\infty} \frac{1}{16^n} \left( \frac{4}{8n+1} - \frac{2}{8n+4} - \frac{1}{8n+5} - \frac{1}{8n+6} \right)[/math]

 

Integer sum:

[math]n \; (0,0,1,1,2,3,3,4,5)[/math]

 

Edited July 10, 2012 by Orion1

[Orion1]

Madhava-Leibniz series: (14th century)

 

[Math]\pi = \sqrt{12} \sum_{k=0}^{\infty} \frac{(-3)^{-k}}{2k + 1}[/Math]

 

Integer sum:

[math]k \; (0,2,4,5,7,8,11,13,14)[/math]

Reference:

Madhava-Leibniz series - Wikipedia

 

Chudnovsky algorithm: (the fastest algorithm in its class)

 

[math]\frac{1}{\pi} = 12 \sum^\infty_{k=0} \frac{(-1)^k (6k)! (13591409 + 545140134k)}{(3k)!(k!)^3 640320^{3k + 3/2}}[/math]

 

Integer sum:

[math]k \; (0,0,0,0,0,0,0,0,0)[/math]

 

[math]n_{19} = 7[/math] at [math]k = 0[/math]

[math]n_{28} = 3[/math] at [math]k = 1[/math]

Mathematica source code:

N[1/(12*Sum[((-1)^k*Factorial[6*k]*(13591409 + 545140134*k))/(Factorial[3*k]*Factorial[k]^3*640320^(3*k + 3/2)), {k, 0, 1}]), 30]

Reference:

Chudnovsky algorithm - Wikipedia

Edited July 11, 2012 by Orion1

[IsaacAsimov]

If the computer is evaluating pi from scratch, such as the case example for the Atn(x) and partial sum functions, then the computer does not know that the answer is pi, so how can it 'know' to round up to the nearest decimal point? and in the case where the solution is only 3.141 at that point in the sum, then there are no further decimal places to round up from!.

 

 

Excellent!, I used to program on an old Apple IIe myself, the basic languages are very similar. The C64 is really easy to program.

 

Have you tried using a Commodore 64 emulator yet?

 

Check out the WinVICE-2.2-x86 link listed in reference 1, and try out the 'warp' option!.

Reference:

WinVICE-2.2-x86 - a Commodore 64 emulator for Windows

 

I downloaded WinVICE, and I found the Commodore 64 screen, but I couldn't get it to work properly.

 

Besides, why would I need a C64 emulator when I have an actual C64?

 

Thank you for your concern.

[doG]

Mathematica source code:

N[1/(12*Sum[((-1)^k*Factorial[6*k]*(13591409 + 545140134*k))/(Factorial[3*k]*Factorial[k]^3*640320^(3*k + 3/2)), {k, 0, 1}]), 30]

Did you know you can paste Mathematica source code directly into WolframAlpha for evaluation?

[IsaacAsimov]

I downloaded WinVICE, and I found the Commodore 64 screen, but I couldn't get it to work properly.

 

Besides, why would I need a C64 emulator when I have an actual C64?

 

Thank you for your concern.

 

I finally got the C64 emulator to work, after I found out that it uses the same keyboard configuration as an actual C64.

---

Pi 16K8 Formula:

 

Here is a perfect pi formula that I have found.

It converges on pi very quickly, and doesn't need to use factorials.

 

 

Here is the program to go with the formula. I wrote it in Structured BASIC on my C64 computer.

You can use values of k from 0 to 31 on the C64 without getting an error.

I got rid of the unnecessary brackets in line 340.

 

It is a very short program.

 

100 REM PI 16K8

120 REM FROM PI FORMULA USING 16^K AND 8

140 REM

160 CALL MAIN

180 END

200 :

220 PROC MAIN

240 .....L=10

260 .....REM PRINT "A",,"S"

280 .....K=-1

300 .....LOOP

320 ..........K=K+1

340 ..........A=1/(16^K)*(4/(8*K+1)-2/(8*K+4)-1/(8*K+5)-1/(8*K+6))

360 ..........S=S+A

380 ..........REM PRINT A,S

400 .....UNTIL K=L

420 .....PI=S

440 .....PRINT:PRINT "PI =";PI

460 ENDPROC

[doG]

Why don't you try implementing the Gauss–Legendre iterative algorithm. It is notable for being rapidly convergent, with only 25 iterations producing 45 million correct digits of π. That should rapidly bring your Commodore 64 to it's knees

[IsaacAsimov]

Here's a pi formula I found in 'A History of Pi' by Petr Beckmann, top of p. 185:

 

[math]\pi[/math] = 24 arctan(1/8)+8 arctan(1/57)+4 arctan(1/239)

 

Here's the program, typed on a C64 computer, that goes with it:

 

10 REM PI24ATN8 FROM 'A HISTORY OF PI'

20 REM TOP OF PAGE 185

30 PI=24*ATN(1/8)+8*ATN(1/57)+4*ATN(1/239)

40 PRINT

50 PRINT "PI =";PI

[John Cuthber]

[math]\pi[/math] =[math]\pi[/math]

[baxtrom]

The quarter circle arc length method, although being intuitive, is problematic since any numerical integration algorithm will converge slowly due to nasty integrand when approaching singularity at x = 1. However, one way of getting around this is to compute the arc length of a 1/8 unit circle, which should be pi/4. Hence, replacing the upper integration limit 1 -> sqrt(0.5) and doubling the constant in front of the integral, the identity becomes

 

[math]\pi = 4 \int_{0}^{\sqrt{1/2}} \sqrt { 1 + \frac{x^2}{1-x^2} } \mathrm{d}x [/math]

 

Using Simpson's rule with 16 subintervals gives pi ~ 3.14160046787255 (correct digits in bold)

with 512 subintervals and Richardson extrapolation: pi ~ 3.14159265358979

 

In comparison, using the area method (computing the area of a unit circle), i.e.

 

[math]\pi = 4 \int_{0}^{1} \sqrt { 1 + x^2 } \mathrm{d}x [/math]

 

..one arrives at pi ~ 3.14159246129562 after 16384 subintervals + Richardson extrapolation.

 

It looks like the 1/8 arc length method is converging faster than the area method.

Edited July 17, 2012 by baxtrom

[IsaacAsimov]

[math]\pi[/math] =[math]\pi[/math]

 

That is complete nonsense. Anybody else who has something worthwhile to say about post #107, please do so.

[John Cuthber]

That is complete nonsense. Anybody else who has something worthwhile to say about post #107, please do so.

No.

It isn't nonsense.

You just didn't understand it.

Essentially, it's just a shorter version of the programme you gave.

What you seem to have missed is that, in order to calculate arctan you need to calculate an infinite series so you might as well calculate one that gives pi directly.

Calculating 3 series is going to be slower (in general) than just calculating 1 so your method, while mathematically amusing, isn't practical.

 

Also, the series sums used in computers are carefully optimised to run quickly and converge well. For example they are more likely to use Chebychev polynomials or this sort of thing

http://en.wikipedia.org/wiki/CORDIC

than Taylor series.

However those series require a prior knowledge of the various constants.

I strongly suspect that, at least for trig functions, one of those is pi.

So, I think your method needs to know what pi is in order to calculate pi.

 

My method does the same but is quicker.

 

Incidentally, even setting that aside, statements of fact are not nonsense, so your first sentence is wrong anyway.

As for "Anybody else who has something worthwhile to say about post #107, please do so"

Well, since this is a discussion forum, if they had anything to say, they would have done so.

So your second sentence is redundant.

In future you might want to not bother with posts that are half wrong and half pointless.

Edited July 18, 2012 by John Cuthber

[CaptainPanic]

!

Moderator Note

John Cuthber, please be so kind to provide a little more explanation to avoid the situation where your posts are misunderstood. This refers specifically to the post of [math]\pi[/math] =[math]\pi[/math].
IsaacAsimov, please don't be provoked immediately and don't respond so strongly - it can only lead to escalation.

[John Cuthber]

Sorry, I thought it was obvious.

[baxtrom]

Did anyone mention trial and error? There's a fancy name for this type of analysis: Monte carlo algorithms!

Imagine you are throwing a dart onto a 2x2 square (area = 4). You are good enough to always hit the square, but them darts land all over the place. Now, a ratio of those darts will fall inside a unit circle concentric with the square (area = pi). That ratio is the same as the ratio of areas: pi / 4. I assume you have a random number generator in your c64. Then you can write a program which does something like this:

 

randomize

n = 0

N = big number!

 

for i = 1 to N

(creates random coordinate point in x = -1..1, y = -1..1)

x = 2*rand-1

y = 2*rand-1

R = sqrt(x^2 + y^2)

 

(checks if [x,y] is within unit circle, if so -> increase n)

if R <= 1 then

n = n + 1

end

end

 

(n/N is approximately pi/4)

pi_approx = 4*n/N;

 

Disclaimer, code may be wrong, didn't test it. Even if it works, it will be a seriously slow method for getting correct decimals of pi

[John Cuthber]

It would be quicker (and give the same result) if you missed out the sqrt function.

 

Also, an equivalent experiment has been done.

http://en.wikipedia.org/wiki/Buffon's_needle#Estimating_.CF.80

[baxtrom]

It would be quicker (and give the same result) if you missed out the sqrt function.

 

True, but then it would require some additional explanation so I skipped it. Same goes with the random coordinates, "2*rand-1" could be replaced with just "rand" since the ratio of areas of a quarter unit circle to a quarter square is also pi/4.

 

Oh, by the way: it's important that the random numbers are uniformly distributed. Typically the rand functions are UD but some generators produce normally distributed random numbers which would screw up the algorithm and produce pi_approx > pi.

[John Cuthber]

That's OK, if the data are normally distributed you can calculate pi from the area under the probability distribution curve.

[baxtrom]

From the department for cumbersome ways of estimating pi

 

[math]\pi = \Gamma^2(1/2)[/math]

 

..where Gamma is the Gamma function,

[math]\Gamma(z) = \int_0^\infty e^{-t} t^{z-1} \mathrm{d}t[/math]

 

I was just going to suggest to use a Stirling approximation* to estimate the gamma function instead of computing the integral, but of course all those approximations contain pi..

 

 

_{*) Yes, these approximations are accurate only for large arguments, but Gamma(0.5) can be rewritten 2 Gamma(1.5) et c}

[Orion1]

I finally got the C64 emulator to work, after I found out that it uses the same keyboard configuration as an actual C64.

The real advantage of using an emulator is that programmers can still program on a modern operating system and can also backup all their programs to a modern hard drive or CD/DVD ROM and still have access to all their programs and games without fear of a damaged floppy disk failure or an eventual hardware failure.

 

It is possible to format the emulator keyboard so that the keyboard symbols match the modern western style keyboards.

 

Locate the ROM file in the C64 folder named 'win_pos.vkm' and rename it to 'win_pos-old.vkm', then simply copy and paste the code below into Wordpad and save it to the C64 folder as 'win_pos.vkm'.

# VICE keyboard mapping file
#
# Modified by Leif Bloomquist on March 23/2007 to
#   give a proper symbolic mapping on US kbds - finally!!!1
#
# A Keyboard map is read in as patch to the current map.
#
# File format:
# - comment lines start with '#'
# - keyword lines start with '!keyword'
# - normal line has 'keysym/scancode row column shiftflag'
#
# Keywords and their lines are:
# '!CLEAR'   			clear whole table
# '!INCLUDE filename'	read file as mapping file
# '!LSHIFT row col'  	left shift keyboard row/column
# '!RSHIFT row col'  	right shift keyboard row/column
# '!VSHIFT shiftkey' 	virtual shift key (RSHIFT or LSHIFT)
# '!UNDEF keysym'    	remove keysym from table
#
# Shiftflag can have the values:
# 0  	key is not shifted for this keysym/scancode
# 1  	key is shifted for this keysym/scancode
# 2  	left shift
# 4  	right shift
# 8  	key can be shifted or not with this keysym/scancode
# 16 	deshift key for this keysym/scancode
# 32 	another definition for this keysym/scancode follows
#
# Negative row values:
# 'keysym -1 n' joystick #1, direction n
# 'keysym -2 n' joystick #2, direction n
# 'keysym -3 0' first RESTORE key
# 'keysym -3 1' second RESTORE key
# 'keysym -4 0' 40/80 column key
# 'keysym -4 1' CAPS (ASCII/DIN) key
#

!CLEAR
!LSHIFT 1 7
!RSHIFT 6 4
!VSHIFT RSHIFT

#0 -1 -1 0 			/*   		(no key)   		*/
1 7 7 8            	/*      	ESC -> Run/Stop 	*/
2 7 0 8            	/*        	1 -> 1        	*/
3 7 3 40   			/*        	2 -> 2        	*/
3 5 6 16   			/*        	@ -> @        	*/
4 1 0 8            	/*        	3 -> 3        	*/
5 1 3 8            	/*        	4 -> 4        	*/
6 2 0 8            	/*        	5 -> 5        	*/
7 2 3 40   			/*        	6 -> 6        	*/
7 6 6 16   			/*        	^ -> ^        	*/
8 3 0 40   			/*        	7 -> 7        	*/
8 2 3 1            	/*        	& -> &        	*/
9 3 3 40   			/*        	8 -> 8        	*/
9 6 1 16   			/*        	* -> *        	*/
10 4 0 40          	/*        	9 -> 9        	*/
10 3 3 1   			/*        	( -> (        	*/
11 4 3 40          	/*        	0 -> 0        	*/
11 4 0 1   			/*        	) -> )        	*/
12 5 3 8   			/*    	Minus -> Minus    	*/
13 6 5 40          	/*    	Equal -> Equal    	*/
13 5 0 16          	/*        	+ -> +        	*/
14 0 0 8   			/*	Backspace -> Del      	*/
15 7 2 8   			/*      	TAB -> Ctrl 		*/
16 7 6 8   			/*        	Q -> Q        	*/
17 1 1 8   			/*        	W -> W        	*/
18 1 6 8   			/*        	E -> E        	*/
19 2 1 8   			/*        	R -> R        	*/
20 2 6 8   			/*        	T -> T        	*/
21 3 1 8   			/*        	Y -> Y        	*/
22 3 6 8   			/*        	U -> U        	*/
23 4 1 8   			/*        	I -> I        	*/
24 4 6 8   			/*        	O -> O        	*/
25 5 1 8   			/*        	P -> P        	*/
26 5 5 1   			/*        	[ -> [        	*/
27 6 2 1   			/*        	] -> ]        	*/
28 0 1 8   			/*   	Return -> Return   	*/
29 7 5 8   			/*	Left Ctrl -> CBM      	*/
30 1 2 8   			/*        	A -> A        	*/
31 1 5 8   			/*        	S -> S        	*/
32 2 2 8   			/*        	D -> D        	*/
33 2 5 8   			/*        	F -> F        	*/
34 3 2 8   			/*        	G -> G        	*/
35 3 5 8   			/*        	H -> H        	*/
36 4 2 8   			/*        	J -> J        	*/
37 4 5 8   			/*        	K -> K        	*/
38 5 2 8   			/*        	L -> L        	*/
39 6 2 40          	/*        	; -> ;        	*/
39 5 5 16          	/*        	: -> :        	*/
40 3 0 33          	/*        	' -> '        	*/
40 7 3 1   			/*        	' -> '        	*/
41 7 1 40          	/*        	` -> Left Arrow   */
41 6 6 1   			/*        	~ -> Pi   		*/
42 1 7 2   			/*   Left Shift -> Left Shift   */
43 6 0 8   			/*        	\ -> Pound    	*/
44 1 4 8   			/*        	Z -> Z        	*/
45 2 7 8   			/*        	X -> X        	*/
46 2 4 8   			/*        	C -> C        	*/
47 3 7 8   			/*        	V -> V        	*/
48 3 4 8   			/*        	B -> B        	*/
49 4 7 8   			/*        	N -> N        	*/
50 4 4 8   			/*        	M -> M        	*/
51 5 7 8   			/*        	, -> ,        	*/
52 5 4 8   			/*        	. -> .        	*/
53 6 7 8   			/*        	/ -> /        	*/
54 6 4 4   			/*  Right Shift -> Right Shift  */
55 6 1 8   			/*   	Grey * -> *        	*/
#56 -1 -1 0        	/* 	Left Alt -> (no key) 	*/
57 7 4 8   			/*    	Space -> Space    	*/
58 7 7 8   			/*	Caps Lock -> Run/Stop 	*/
59 0 4 8   			/*   		F1 -> F1   		*/
60 0 4 1   			/*   		F2 -> F2   		*/
61 0 5 8   			/*   		F3 -> F3   		*/
62 0 5 1   			/*   		F4 -> F4   		*/
63 0 6 8   			/*   		F5 -> F5   		*/
64 0 6 1   			/*   		F6 -> F6   		*/
65 0 3 8   			/*   		F7 -> F7   		*/
66 0 3 1   			/*   		F8 -> F8   		*/
#67 -1 -1 0        	/*   		F9 -> (no key) 	*/
#68 -1 -1 0        	/*      	F10 -> (no key) 	*/
#69 -1 -1 0        	/* 	Num Lock -> (no key) 	*/
#70 -1 -1 0        	/*  Scroll Lock -> (no key) 	*/
#71 -1 -1 0        	/* 	Numpad 7 -> (no key) 	*/
#72 -1 -1 0        	/* 	Numpad 8 -> (no key) 	*/
#73 -1 -1 0        	/* 	Numpad 9 -> (no key) 	*/
#74 -1 -1 0        	/* 	Numpad - -> (no key) 	*/
#75 -1 -1 0        	/* 	Numpad 4 -> (no key) 	*/
#76 -1 -1 0        	/* 	Numpad 5 -> (no key) 	*/
#77 -1 -1 0        	/* 	Numpad 6 -> (no key) 	*/
#78 -1 -1 0        	/* 	Numpad + -> (no key) 	*/
#79 -1 -1 0        	/* 	Numpad 1 -> (no key) 	*/
#80 -1 -1 0        	/* 	Numpad 2 -> (no key) 	*/
#81 -1 -1 0        	/* 	Numpad 3 -> (no key) 	*/
#82 -1 -1 0        	/* 	Numpad 0 -> (no key) 	*/
#83 -1 -1 0        	/* 	Numpad . -> (no key) 	*/
#84 -1 -1 0        	/*   	SysReq -> (no key) 	*/
#85 -1 -1 0        	/*   		85 -> (no key) 	*/
#86 -1 -1 0        	/*   		86 -> (no key) 	*/
#87 -1 -1 0        	/*      	F11 -> (no key) 	*/
#88 -1 -1 0        	/*      	F12 -> (no key) 	*/
89 6 3 8   			/* 		Home -> CLR/HOME 	*/
90 0 7 1   			/*   		Up -> CRSR UP  	*/
#91 -1 -1 0        	/* 		PgUp -> (no key) 	*/
92 0 2 1   			/* 		Left -> CRSR LEFT	*/
93 0 2 8   			/*    	Right -> CRSR RIGHT   */
#94 -1 -1 0        	/*      	End -> (no key) 	*/
95 0 7 8   			/* 		Down -> CRSR DOWN	*/
#96 -1 -1 0        	/*   	PgDown -> (no key) 	*/
97 0 0 1   			/*      	Ins -> Shift-Del (Ins)*/
98 0 0 8   			/*      	Del -> Del      	*/
#99 -1 -1 0        	/* Numpad Enter -> (no key) 	*/
#100 -1 -1 0   		/*   Right Ctrl -> (no key) 	*/
#101 -1 -1 0   		/*    	Pause -> (no key) 	*/
#102 -1 -1 0   		/*   	PrtScr -> (no key) 	*/
#103 -1 -1 0   		/* 	Numpad / -> (no key) 	*/
#104 -1 -1 0   		/*	Right Alt -> (no key) 	*/
#105 -1 -1 0   		/*    	Break -> (no key) 	*/
106 7 5 8          	/*   Left Win95 -> CBM      	*/
#107 -1 -1 0   		/*  Right Win95 -> (no key) 	*/

#
# Joystick 1
#
#KP_0 -1 0
#KP_1 -1 1
#KP_2 -1 2
#KP_3 -1 3
#KP_4 -1 4
#KP_5 -1 5
#KP_6 -1 6
#KP_7 -1 7
#KP_8 -1 8
#KP_9 -1 9

#
# Joystick 2
#
#w -2 7
#e -2 8
#r -2 9
#s -2 4
#d -2 5
#f -2 6
#x -2 1
#c -2 2
#v -2 3
#space -2 0

# Restore key mappings
91 -3 0 		/* 		PgUp -> RESTORE  	*/

Note that the Commodore 64 can only display eight decimal places and the basic commands can only see eight decimal places. The C64 is capable of analyzing nine decimal places, however at ten decimal places, an infinity is produced and no further decimal evaluation is possible.

 

So for equation iterations that produce more decimal places than nine, are incapable of being further evaluated by the Commodore 64.

 

Note that ^ is code for [math]\pi[/math] and ~ is code for ^ under the VICE emulator copy/paste function.

10 clr
20 n1 = 12
30 v1 = val(mid$(str$(^),n1,1))
40 print v1:print ^
50 pi = ^
60 n2 = 0:e1 = -9
70 v1 = val(mid$(str$(pi),11,1))
80 if v1 = 7 then goto 110
90 pi = pi + 10~(e1):n2 = n2 + 1
100 goto 70
110 print n2""pi
120 end

Reference:

The VICE Emulator

 

Edited July 19, 2012 by Orion1

[IsaacAsimov]

Thank you for the information, Orion1. I might need it someday.

[IsaacAsimov]

Note that the Commodore 64 can only display eight decimal places and the basic commands can only see eight decimal places. The C64 is capable of analyzing nine decimal places, however at ten decimal places, an infinity is produced and no further decimal evaluation is possible.

 

So for equation iterations that produce more decimal places than nine, are incapable of being further evaluated by the Commodore 64.

 

Note that ^ is code for and ~ is code for ^ under the VICE emulator copy/paste function.

10 clr

20 n1 = 12

30 v1 = val(mid$(str$(^),n1,1))

40 print v1:print ^

50 pi = ^

60 n2 = 0:e1 = -9

70 v1 = val(mid$(str$(pi),11,1))

80 if v1 = 7 then goto 110

90 pi = pi + 10~(e1):n2 = n2 + 1

100 goto 70

110 print n2""pi

120 end

 

What is this program supposed to accomplish?

[Orion1]

What is this program supposed to accomplish?[/Quote]

It is a test of the computational numerical limits of a Commodore 64 computer.

 

For equation iterations that produce more decimal places than nine, are incapable of being further evaluated by the Commodore 64.

[IsaacAsimov]

Using My Calculator:

 

Back To The Future: 88 mi/h x 1.609 km/mi = 141.592 km/h [math]\simeq[/math] 142 km/h

 

Pi'th root of pi: [math]\sqrt[\pi]{pi}=1.43962 [/math]

 

Pi to the power of pi: [math] \pi^\pi = 36.46216 [/math]