Computing PI
Author
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.
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]
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]
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.
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; }}
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.
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.
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.
http://fabrice.bellard.free.fr/pi/index.html
algorithm that gets n digits in O(n^2):
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;} 
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
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
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)
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!
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!
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
