Methods of Estimating Pi

[IsaacAsimov]

How about (2/3)(4/5)(6/7)..... instead of doing the numerator and denominator separately?

 

By golly, that's a great idea! Thanks!

[khaled]

[math]f(m) = 2 \times 4 \times 6 \times 8 \times ... \times (2m-2)[/math]

 

[math]= 2 \times [ 1 \times 2 \times 3 \times 4 \times ... \times (m-1) ][/math]

 

[math]= 2 \times (m-1)![/math]

Edited May 23, 2012 by khaled

[imatfaal]

[math]f(m) = 2 \times 4 \times 6 \times 8 \times ... \times (2m-2)[/math]

 

[math]= 2 \times [ 1 \times 2 \times 3 \times 4 \times ... \times (m-1) ][/math]

 

[math]= 2 \times (m-1)![/math]

 

Nope again. Try putting some figures in before you post

 

[math]= (2 \times 4 \times 6)=48[/math]

[math]= 2 \times (1 \times 2 \times 3)=12[/math]

 

 

You are dividing the contents of the bracket by 2 multiple times - ie for each value inside; so you need to multiply by 2 as many times

 

[math](2 \times 4 \times 6)=2 \times (1 \times 4 \times 6)=2 \times 2 \times (1 \times 2 \times 6)=2^3 \times (1 \times 2 \times 3)[/math]

 

[math] f(m) = 2 \times 4 \times 6 \times 8 \times ... \times (2m-2) =2^m \times (m-1)! [/math]

[hypervalent_iodine]

!

Moderator Note

As per my previous modnote, this thread is now the culmination of a number of threads written by IsaacAsimov on different methods for approximating pi. Although normally we would close threads considered as duplicates, it's obvious that Isaac has put in a lot of work in his calculations and so I thought it would be better if I instead merged them all into an all-in-one SUPER pi thread.

Isaac, as mentioned, I would like it if you could please post anything further related to approximations of pi into this thread rather than posting more and more new ones.

[Xittenn]

I'm a little confused why the title includes 'sneaky' . .. . it sort of implies that a more straight forward method of compiling an accurate answer exists, which which it doesn't. I believe simply "Methods of Estimating Pi" would have been most appropriate, however less fun . . . .

 

or "Numerical Methods of Estimating Pi"

Edited May 25, 2012 by Xittenn

[hypervalent_iodine]

For the record, I didn't name the thread.

[IsaacAsimov]

Here's a program I wrote in Structured BASIC on my C64 computer.

 

100 REM COMPUTING PI USING TRAPEZOIDS

120 CALL INIT

140 CALL MAIN

160 CALL OUTPUT

180 END

20O :

220 PROC INIT

240 .....TI$="000000"

260 .....A=0

280 .....B=1

300 .....X%=1

320 .....N=10

340 ENDPROC

360 :

380 PROC MAIN

400 .....DX=(B-A)/N

420 .....PRINT:PRINT "DX =";DX

440 .....X=-DX

460 .....I=-1

480 .....PRINT:PRINT " I ";" X ";" Y ";" T ":PRINT

500 .....LOOP

520 ..........I=I+1

540 ..........X=X+DX

560 ..........D=1-X*X

580 ..........REM

600 ..........IF D<0

620 ...............D=0

640 ..........ENDIF

660 ..........REM

680 ..........Y=2*SQR(D)

700 ...........IF I=0 OR I=N

720 ................Y=Y/2

740 ...........ENDIF

760 ...........T=T+Y

780 ...........PRINT I;X;Y;T

800 .....UNTIL I=N

820 .....PI=DX/2*T*4

840 ENDPROC

860 :

880 PROC OUTPUT

900 .....PRINT:PRINT "PI =";PI

920 .....PRINT:PRINT "TIME TAKEN = ";TI$

940 ENDPROC

[IsaacAsimov]

 

Here's a Structured BASIC program I wrote on my C64 computer that computes the even/odd pi formula shown above.

 

100 REM EVEN/ODD PI FORMULA

110 REM DOTS REPRESENT INDENTATIONS

120 CALL INIT

140 CALL MAIN

160 CALL OUTPUT

180 END

190 :

200 PROC INIT

220 .....TI$="000000"

240 .....EV=0:OD=1

260 .....P=1

280 .....L=500

300 .....EP=2*L-2

320 ENDPROC

340 :

360 PROC MAIN

380 .....M=1

400 .....REM PRINT"EV";" OD ";"F",,"P"

420 .....PRINT

440 .....LOOP

460 ..........M=M+1

480 ..........EV=2*M-2

500 ..........OD=2*M-1

520 ..........F=EV/OD

540 ..........P=P*F

560 ..REM PRINT EV;OD;F;P

580 .....UNTIL EV=EP

600 .....REM

620 .....A=P*SQR(2*M)

640 .....PI=2*A*A

660 .....REM

680 ENDPROC

700 :

720 PROC OUTPUT

740 .....PRINT:PRINT"A =";A

760 .....PRINT:PRINT"PI =";PI

780 .....PRINT:PRINT"TIME TAKEN = ";TI$

800 ENDPROC

[ewmon]

I don't see the following formula among the posts here:

 

π ~ n∙sin(360/n)/2

as n → ∞ (start with 4 and go by factors of 2).

 

This method uses regular polygons inscribed within a circle. It starts with an inscribed square, then an octagon, 16-gon, 32-gon, etc. It's interesting to note that each new n reduces the error by about 25%. This means that each reiteration adds to the approximation about ¾ of the segment between the chord and the arc.

 

[IsaacAsimov]

There's a formula very similar to the formula you have given in post #21 of this thread. Check it out.

[IsaacAsimov]

By Arclength Approximation Riemann Form

 

[math] f(x) = \sqrt { 1 - x^2 }\, ; \; f'(x) = - \frac{x}{\sqrt{1-x^2}}\, ; \; f''(x) = - \frac{1}{\sqrt{(1-x^2)^3}}\, ; \; f'''(x) = - \frac{3x}{\sqrt{(1-x^2)^5}}[/math]

 

 

 

[math] S = \frac {\pi}{2} = \int_{0}^{1} \sqrt { 1 + [f'(x)]^2 } \, dx = \int_{0}^{1} \sqrt { 1 + [- \frac{x}{\sqrt{1-x^2}}]^2 } \, dx = \int_{0}^{1} \sqrt { 1 + \frac{x^2}{1-x^2} } \, dx [/math]

 

 

 

[math] \Delta x = \frac{b - a}{n} = \frac{1}{n} [/math]

 

[math] x_i^* = a + i \cdot \Delta x = \frac{i}{n} [/math]

 

[math] x_{i-1}^* = a + (i - 1) \cdot \Delta x = \frac{i - 1}{n} [/math]

 

[math] \bar{x} = \frac{x_i^* + x_{i-1}^*}{2} = \frac{(\frac{i}{n}) + (\frac{i - 1}{n})}{2} = \frac{2i - 1}{2n} [/math]

 

 

 

[math] M_n \approx \sum_{i=1}^n f( \bar{x} ) \cdot \Delta x = \sum_{i=1}^n \sqrt { 1 + \frac{(\bar{x})^2}{1-(\bar{x})^2} } \cdot \frac{1}{n} = \frac{1}{n} \sum_{i=1}^n \sqrt { 1 + \frac{(\frac{2i - 1}{2n})^2}{1-(\frac{2i - 1}{2n})^2} } = \frac{1}{n} \sum_{i=1}^n \sqrt{\frac{4n^2}{4n^2 - 4i^2 + 4i - 1 }}[/math]

 

 

 

[math] error = |S - M_n| \, ; \; K_2 = max | f''(x) | \, ; \; 0 \leq n = \sqrt{\frac{K_2(b - a)^3}{24(|S - M_n|)}} [/math]

 

It looks like you put a lot of effort into this post. The formulas you use are very complicated.

 

Unfortunately, I don't know what any of it means, except it has something to do with pi, which is worthy in itself.

[Xittenn]

As always I made a slight oversight. Not that anything is technically wrong with what was posted, but these three derivatives are useless:

 

[math] f'(x) = - \frac{x}{\sqrt{1-x^2}}\, ; \; f''(x) = - \frac{1}{\sqrt{(1-x^2)^3}}\, ; \; f'''(x) = - \frac{3x}{\sqrt{(1-x^2)^5}} [/math]

 

What I should have done is given the derivatives of the function for which we are integrating:

 

[math] g(x) = \sqrt { 1 + [f'(x)]^2 } [/math]

 

where the importance of these are in calculating the error of the approximation using the equations on the final line.

 

At any rate I'll write up an explanation tomorrow, when I'm at home studying. I hope you are somewhat familiar with calculus? It's really not too complicated. The important part is:

 

[math] M_n \approx \frac{1}{n} \sum_{i=1}^n \sqrt{\frac{4n^2}{4n^2 - 4i^2 + 4i - 1 }} [/math]

Edited June 1, 2012 by Xittenn

[baxtrom]

It looks like you put a lot of effort into this post. The formulas you use are very complicated.

 

Unfortunately, I don't know what any of it means, except it has something to do with pi, which is worthy in itself.

 

It's an algorithm for approximating pi using numerical integration (midpoint rule) to obtain the arc length of a quarter unit circle (which equals pi/2). For a given function the arc length can be computed using integral calculus. In this case, however, the integrand is not very well behaved close to x = 1 so accuracy may suffer. For a more practical method using integral calculus I would suggest instead computing the area of a quarter unit circle = pi/4.

[ewmon]

Interesting. I made a polygon perimeter method to compare to my polygon area method in post #59, and it's estimations using n-gons equal the area estimations using 2n-gons.

Edited June 1, 2012 by ewmon

[baxtrom]

I don't see the following formula among the posts here:

 

π ~ n∙sin(360/n)/2

as n → ∞ (start with 4 and go by factors of 2).

 

 

I guess a problem with that formula is that it requires implicit knowledge of pi beforehand. This is more clear if you rewrite the limit in terms of radians instead. Your formula would then become

 

[math]\pi = \lim_{n \to \infty} \frac{n \sin(\frac{2 \pi}{n})}{2} = \text{(sine expansion for small arguments)} = \lim_{n \to \infty} \frac{\frac{2 n \pi}{n}}{2} \ldots[/math]

[Xittenn]

In this case, however, the integrand is not very well behaved close to x = 1 so accuracy may suffer. For a more practical method using integral calculus I would suggest instead computing the area of a quarter unit circle = pi/4.

 

Which is definitely true of the original improper integral and would be true of the [math] M_n [/math] approximation if [math] i = n + \frac{1}{2} [/math] produced a whole number given whole numbered n--do correct me if I'm wrong! The original thread title under which I posted required a quarter circle [math] \frac{2pi}{4} [/math].

 

** also where [math] i \not > n [/math]

Edited June 1, 2012 by Xittenn

[Xittenn]

So yeah, explanation . .. . . . .

 

 

First we realize that the following function [math] f(x) [/math] is simply the equation for the unit circle, and we do in fact require the derivative [math] f'(x) [/math] to plug into the subsequent integral equation:

 

[math] f(x) = \sqrt { 1 - x^2 }\, ; \; f'(x) = - \frac{x}{\sqrt{1-x^2}} [/math]

 

 

Now by definition [math] 2\pi [/math] is the length of the arc that is the circumference of the unit circle, and so [math] \frac{\pi}{2} [/math] is [math] \frac{1}{4} [/math] of the circumference of the unit circle. We use only a fraction of the unit circle to simplify the overall problem. According to a simple theorem we can derive the arclength of a curve by breaking it down into an infinite number of pieces such that if we add the length of each 'chord' we will find the length of the arc. We can say that a small change in [math] x [/math] will create a small change in [math] y = f(x) [/math] such that the length of the chord can be written as follows:

 

[math] |P_{i-1}P_{i}| = \sqrt{(\Delta x)^2 + (\Delta y)^2 } = \sqrt{[\Delta x]^2 +[f'(x^*_i) \Delta x]^2} = \sqrt{1+[f'(x^*_i)]^2} \sqrt{[\Delta x]^2} = \sqrt{1+[f'(x^*_i)]^2} \Delta x [/math]

 

which by taking the limit of the infinite number of pieces of the arc, or by adding the chords as they become infinitely small and taking their limit:

 

[math] L = \lim_{n \to \infty} \sum^n_{i=1} \sqrt{1+[f'(x^*_i)]^2} \Delta x = \int^b_a \sqrt{1+[f'(x)]^2} dx = S[/math]

 

This is the arclength formula and is what I'm using for the circumference of the quarter circle. I'm not trying to be a smart ass, I'm just trying to give a proper elaboration that might help you see the correlation between the weird stuff and what you are already doing, without posting something that is inaccurate.

 

 

So moving on we simply plug in [math] f'(x) [/math] from above into our arclength integral. If you really wanted to you could solve this integral using limits . . . . 'if' !

 

[math] S = \frac {\pi}{2} = \int_{0}^{1} \sqrt { 1 + [f'(x)]^2 } \, dx = \int_{0}^{1} \sqrt { 1 + [- \frac{x}{\sqrt{1-x^2}}]^2 } \, dx = \int_{0}^{1} \sqrt { 1 + \frac{x^2}{1-x^2} } \, dx [/math]

 

 

But we don't want to do this. We want to reinterpret this integral as a Riemann Sum so that you can use your C64 to solve for pi in basic, because this seems to be what you love to do??

 

The basic form for taking an integral and reinterpreting it as a Riemann Sum is as follows:

 

[math] \int^b_a f(x) dx [/math]

 

We make the following analogies between the two interpretations of our structure:

 

width of trapezoid such that [math] n [/math] is the number of divisions or 'cuts':

 

[math] dx = \Delta x = \frac{b - a}{n} = \frac{1}{n} [/math]

 

 

left trapezoid edge as it meets the curve for mid approximation:

 

[math] x = x_i^* = a + i \cdot \Delta x = \frac{i}{n} [/math]

 

 

right trapezoid edge as it meets the curve for mid approximation:

 

[math] x_{i-1}^* = a + (i - 1) \cdot \Delta x = \frac{i - 1}{n} [/math]

 

 

the average of the left and right trapezoids as the midpoint meets the curve for mid approximation:

 

[math] \bar{x} = \frac{x_i^* + x_{i-1}^*}{2} = \frac{(\frac{i}{n}) + (\frac{i - 1}{n})}{2} = \frac{2i - 1}{2n} [/math]

 

 

the integral is replaced by the summation symbol and as it is a mid approximation i=1 and n=n:

 

[math] \int = \sum^n_{i=1} [/math]

 

 

We then take [math] \bar{x} [/math] and plug it into the summation via [math] f(x^*_i) [/math] and substitute in [math] \Delta x. [/math] [math] \frac{1}{n} [/math] moves outside of the sum by summation rules and so on . . . solve for the approximation.

 

[math] M_n \approx \sum_{i=1}^n f( \bar{x} ) \cdot \Delta x = \sum_{i=1}^n \sqrt { 1 + \frac{(\bar{x})^2}{1-(\bar{x})^2} } \cdot \frac{1}{n} = \frac{1}{n} \sum_{i=1}^n \sqrt { 1 + \frac{(\frac{2i - 1}{2n})^2}{1-(\frac{2i - 1}{2n})^2} } = \frac{1}{n} \sum_{i=1}^n \sqrt{\frac{4n^2}{4n^2 - 4i^2 + 4i - 1 }} [/math]

 

 

The stuff on the bottom is to find the error bounds with given n, but as mentioned before it requires the first, second, and third derivatives of the function found within the integral or arbitrarily [math] g(x) [/math]. I'll gladly post these or the [math] K_2 [/math] if you want. The important part though is the summation or the last variant thereof:

 

[math] \frac{1}{n} \sum_{i=1}^n \sqrt{\frac{4n^2}{4n^2 - 4i^2 + 4i - 1 }} [/math]

 

which if you make a little program and make n sufficiently large you will see [math] \frac{pi}{2} [/math]. Simple enough to double this value et la voila, c'est tout! If I made any mistakes, sorry, I'm a noob and it's in my nature . . . .

[IsaacAsimov]

The important part though is the summation or the last variant thereof:

 

 

which if you make a little program and make n sufficiently large you will see

 

Here is a program I wrote in Structured BASIC on my C64 computer which computes the above formula.

 

 

 

100 REM PI USING 4

120 :

140 CALL INIT

160 CALL MAIN

180 CALL OUTPUT

190 END

200 :

220 PROC INIT

240 .....N=10

300 ENDPROC

320 :

340 PROC MAIN

350 .....A=1/N

352 .....PRINT " I ";" B ",," S "

360 .....I=0

380 .....LOOP

400 ..........I=I+1

420 ..........B=SQR((4*N*N)/(4*N*N-4*I*I+4*I-1))

440 ..........S=S+B

460 ..........PRINT I;B;S

480 .....UNTIL I=N

490 .....P1=A*S

500 .....PI=P1*2

520 ENDPROC

540 :

560 PROC OUTPUT

580 .....PRINT:PRINT "PI =";PI

60O ENDPROC

[IsaacAsimov]

In 1910, the Indian mathematician Srinivasa Ramanujan found several rapidly converging infinite series of pi, including

 

 

which computes a further eight decimal places of pi with each term in the series. His series are now the basis for the fastest algorithms currently used to calculate pi.

 

 

 

Here is a program I wrote in Structured BASIC with my C64 computer which computes the above formula:

 

100 REM PI RAMA

110 REM FROM SRINIVASA RAMANUJAN

120 REM

140 CALL INIT

160 CALL MAIN

180 CALL OUTPUT

200 END

220 :

240 PROC INIT

260 .....L=3

280 ENDPROC

300 :

320 PROC MAIN

340 .....A=2*SQR(2)/9801

360 .....K=-1

380 .....LOOP

400 ..........K=K+1

420 ..........P=4*K

440 ..........CALL FACT

460 ..........F1=F

480 ..........P=K

500 ..........CALL FACT

520 ..........F2=F

540 ..........B=F1*(1103+26390*K)/(F2^4*396^(4*K))

560 ..........S=S+B

580 ..........REM PRINT S

600 .....UNTIL K=L

620 .....P1=A*S

640 .....PI=1/P1

660 ENDPROC

680 :

700 PROC OUTPUT

720 .....PRINT:PRINT "PI =";PI

740 ENDPROC

760 :

780 PROC FACT

800 .....I=P

820 .....F=1

840 .....LOOP

860 ..........F=F*I

880 ..........I=I-1

900 ..........REM PRINT I;F

920 .....UNTIL I<=0

940 .....IF I=-1

960 ..........F=1

980 .....ENDIF

990 ENDPROC

[IsaacAsimov]

In 1989, the Chudnovsky brothers correctly computed pi to over a 1 billion decimal places on the supercomputer IBM 3090 using the following variation of Ramanujan's infinite series of pi:

 

 

Here is a computer program I wrote in Structured BASIC on my C64 computer which computes the above formula (includes a procedure for finding factorials).

 

100 REM PI CHUD

110 REM FROM THE CHUDNOVSKY BROTHERS

112 REM DOTS REPRESENT INDENTATIONS

120 REM

140 CALL INIT

160 CALL MAIN

180 CALL OUTPUT

200 END

220 :

240 PROC INIT

260 .....L=1

280 ENDPROC

300 :

320 PROC MAIN

340 .....A=12

360 .....K=-1

380 .....LOOP

400 ..........K=K+1

420 ..........P=6*K

440 ..........CALL FACT

460 ..........F1=F

480 ..........P=3*K

500 ..........CALL FACT

520 ..........F2=F

530 ..........P=K

532 ..........CALL FACT

534 ..........F3=F

540 ..........B=(-1)^K*F1*(13591409+545140134*K)/(F2*F3^3*640320^(3*K+3/2))

560 ..........S=S+B

580 ..........REM PRINT S

600 .....UNTIL K=L

620 .....P1=A*S

640 .....PI=1/P1

660 ENDPROC

680 :

700 PROC OUTPUT

720 .....PRINT:PRINT "PI =";PI

740 ENDPROC

760 :

780 PROC FACT

800 .....I=P

820 .....F=1

840 .....LOOP

860 ..........F=F*I

880 ..........I=I-1

900 ..........REM PRINT I;F

920 .....UNTIL I<=0

940 .....IF I=-1

960 ..........F=1

980 .....ENDIF

990 ENDPROC

[IsaacAsimov]

Pi 16k Formula:

 

Here is the most perfect pi formula I have ever found.

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

 

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

 

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.

It is a very short program.

 

100 REM PI 16K

120 REM FROM PI FORMULA USING 16^K

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

[Orion1]

This simple trigonometry equation trumps all the hyperbole calculus equations stated on this thread so far.

 

Calculation for pi:

[math]4 \cdot \arctan(1) = \pi[/math]

Edited July 7, 2012 by Orion1

[mathematic]

This simple trigonometry equation trumps all the hyperbole calculus equations stated on this thread so far.

 

Calculation for pi:

[math]4 \cdot \arctan(1) = \pi[/math]

 

How would you calculate arctan(1)?

[Orion1]

How would you calculate arctan(1)?[/Quote]

[math]\arctan(z) = \sum_{n=0}^\infty \frac {(-1)^n z^{2n+1}} {2n+1} \; \qquad | z | \le 1 \qquad z \neq i,-i[/math]

 

[math]4 \arctan(1) = 4 \cdot \sum_{n=0}^\infty \frac {(-1)^n 1^{2n+1}} {2n+1} = \pi[/math]

Reference:

Inverse trigonometric functions - infinite series - Wikipedia

[Bignose]

[math]\arctan(z) = \sum_{n=0}^\infty \frac {(-1)^n z^{2n+1}} {2n+1} \; \qquad | z | \le 1 \qquad z \neq i,-i[/math]

 

[math]4 \arctan(1) = 4 \cdot \sum_{n=0}^\infty \frac {(-1)^n 1^{2n+1}} {2n+1} = \pi[/math]

Reference:

Inverse trigonometric functions - infinite series - Wikipedia

 

Ok, but that Ramanujan series converges much faster than this one. That's kind of the point of this thread.