Sign in to follow this  
simon10k

Computing PI

Recommended Posts

I know PI is approx 3.142... but how did "they" work that out? I don't think it was with a string and a ruler, unless it was a very big circle. I heard somewhere that there was a computer program working out the defininate value of PI, that it is a very long number and doesn't seem to have an end. How could I conjure up something in C++ to compute the value of PI without a constant, or just typing it -for learning purposes.

Share this post


Link to post
Share on other sites
Possibly the easiest way to compute PI is using power series expansions. There are more efficient methods, but they tend to be quite a bit more complex. A good overview can be found on Wikipedia:
http://en.wikipedia.org/wiki/History_of_numerical_approximations_of_%CF%80

If you are familiar with Taylor series, you can use the expansion of arctan(x) with Machin's formula (PI/4 = 4*atan(1/5) - atan(1/239)).

[Edited by - Daniel Talagrand on May 24, 2006 5:36:04 AM]

Share this post


Link to post
Share on other sites
A popular approximation - one which is precisely as accurate as you want it to be - is:

4/1 - 4/3 + 4/5 - 4/7 + 4/9 ...

and so on.

I hacked and tested a function to do this. Of course, it can be improved upon. If you want to try to implement it yourself, then ignore the code below.


float pi(float precision)
{
float accumulator = 0, numerator = 4, denominator = 1;
while (1)
{
float term = numerator / denominator;
accumulator += term;
if (abs(term) < precision)
return accumulator;
numerator *= -1;
denominator += 2;
}
}

Share this post


Link to post
Share on other sites
That approximation is even simpler than the one I mentioned, but there is a significant issue with it: the rate of convergence.

Using Machin's formula, if I'm not mistaken, you get around an extra digit of PI per _iteration_. On the other hand, you'd need to do millions upon millions of iterations to obtain anything reasonably accurate using the 4/1-4/3+... series.

Both algorithms are very similar in that they are derived using the Taylor series for atan(x). In this case, PI = 4*atan(1) is used.

Share this post


Link to post
Share on other sites
You can find quite a few different series for pi on this page. The series from Ramanujan are particulary mysterious and have fast rates of convergence.

Share this post


Link to post
Share on other sites
http://fabrice.bellard.free.fr/pi/index.html

algorithm that gets n digits in O(n^2):

/*
* Computation of the n'th decimal digit of \pi with very little memory.
* Written by Fabrice Bellard on January 8, 1997.
*
* We use a slightly modified version of the method described by Simon
* Plouffe in "On the Computation of the n'th decimal digit of various
* transcendental numbers" (November 1996). We have modified the algorithm
* to get a running time of O(n^2) instead of O(n^3log(n)^3).
*
* This program uses mostly integer arithmetic. It may be slow on some
* hardwares where integer multiplications and divisons must be done
* by software. We have supposed that 'int' has a size of 32 bits. If
* your compiler supports 'long long' integers of 64 bits, you may use
* the integer version of 'mul_mod' (see HAS_LONG_LONG).
*/


#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/* uncomment the following line to use 'long long' integers */
/* #define HAS_LONG_LONG */

#ifdef HAS_LONG_LONG
#define mul_mod(a,b,m) (( (long long) (a) * (long long) (b) ) % (m))
#else
#define mul_mod(a,b,m) fmod( (double) a * (double) b, m)
#endif

/* return the inverse of x mod y */
int inv_mod(int x,int y) {
int q,u,v,a,c,t;

u=x;
v=y;
c=1;
a=0;
do {
q=v/u;

t=c;
c=a-q*c;
a=t;

t=u;
u=v-q*u;
v=t;
} while (u!=0);
a=a%y;
if (a<0) a=y+a;
return a;
}

/* return (a^b) mod m */
int pow_mod(int a,int b,int m)
{
int r,aa;

r=1;
aa=a;
while (1) {
if (b&1) r=mul_mod(r,aa,m);
b=b>>1;
if (b == 0) break;
aa=mul_mod(aa,aa,m);
}
return r;
}

/* return true if n is prime */
int is_prime(int n)
{
int r,i;
if ((n % 2) == 0) return 0;

r=(int)(sqrt(n));
for(i=3;i<=r;i+=2) if ((n % i) == 0) return 0;
return 1;
}

/* return the prime number immediatly after n */
int next_prime(int n)
{
do {
n++;
} while (!is_prime(n));
return n;
}

int main(int argc,char *argv[])
{
int av,a,vmax,N,n,num,den,k,kq,kq2,t,v,s,i;
double sum;

if (argc<2 || (n=atoi(argv[1])) <= 0) {
printf("This program computes the n'th decimal digit of \\pi\n"
"usage: pi n , where n is the digit you want\n"
);
exit(1);
}

N=(int)((n+20)*log(10)/log(2));

sum=0;

for(a=3;a<=(2*N);a=next_prime(a)) {

vmax=(int)(log(2*N)/log(a));
av=1;
for(i=0;i<vmax;i++) av=av*a;

s=0;
num=1;
den=1;
v=0;
kq=1;
kq2=1;

for(k=1;k<=N;k++) {

t=k;
if (kq >= a) {
do {
t=t/a;
v--;
} while ((t % a) == 0);
kq=0;
}
kq++;
num=mul_mod(num,t,av);

t=(2*k-1);
if (kq2 >= a) {
if (kq2 == a) {
do {
t=t/a;
v++;
} while ((t % a) == 0);
}
kq2-=a;
}
den=mul_mod(den,t,av);
kq2+=2;

if (v > 0) {
t=inv_mod(den,av);
t=mul_mod(t,num,av);
t=mul_mod(t,k,av);
for(i=v;i<vmax;i++) t=mul_mod(t,a,av);
s+=t;
if (s>=av) s-=av;
}

}

t=pow_mod(10,n-1,av);
s=mul_mod(s,t,av);
sum=fmod(sum+(double) s/ (double) av,1.0);
}
printf("Decimal digits of pi at position %d: %09d\n",n,(int)(sum*1e9));
return 0;
}


Share this post


Link to post
Share on other sites
There is also the "random dart throw" or "random sticks" version of estimating pi.
They are very very very slow to converge, but are very simple to construct.

You inscribe a circle in a square, and "throw darts" randomly distributed into the square,
then check the ratio of ones outside the circle, to the ones inside inorder to get PI.
since the square is (2r)^2 area and the circle is (pi*r^2) area,
pi / 4 = # darts in the circle / #darts

Share this post


Link to post
Share on other sites
Another "random" method of computing pi is to generate two large random numbers and check if they are relatively prime - and repeat many many times. The percentage of those pairs that are relatively prime is about 6/(pi^2)

Share this post


Link to post
Share on other sites
Think about it for a bit. ;-) There are lots of good, established algorithms, but it might be fun to come up with your own.

The obvious method to me back in the day was to find the area of a unit circle. This is an integral which can be approximated by a sum. In rectangular coordinates, each term of the sum involves a square root, which I solved using what I later found out was Newton's method. Pretty straightforward.

The earliest methods of the Greeks, as far as I know, involved circumscribed n-gons. You could come up with a recursive formula for that which you could run to some depth of recursion.

Go make something up!

Share this post


Link to post
Share on other sites
Try this!
#include <stdio.h>

int a=10000,b,c=35000,d,e,f[35001],g;
int main(){
for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),fflush(stdout),e=d%a)
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}



Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this