Printing Pi up to 10,000 digits using C-programming
This program, written in C, allows to to print π to however many decimal places you want and also tells you at the end the computing time.
The original source code was gleaned from the following link (piclassic.c), and was written by Pascal Sebah.
http://numbers.computation.free.fr/Constants/constants.html
This algorithm utilizes the extremely efficient and rapidly converging Machin-like method and formulae to find the value of Pi(π) accurate to up to 10000 decimal places.
When running the program in newer versions of C, however, I found that despite the program working correctly, it printed pi to 10000 decimal places arbitrarily and also seemed to show some runtime errors. I have fixed these, rearranged the program from a single source file to two plus a header to make it neater, changed the function layouts (did old versions of C call them differently?), and also altered it to allow the user to choose the number of digits they would like pi to be printed to.
Original base code however, goes to Pascal Sebah.
If you are using this for an assignment in school/university/etc. then credit both me and Pascal Sebah and for goodness sake change some of the coding or the layout of it so that it's at least some of your own work. At the very least, alter the names of the project, main, functions and header source files. And oh, if it doesn't work, it's because you didn't alter the name of the header file, nor the coding to reflect that.
Remember: Crediting is your friend. I'm not trying to be nasty here, it really is for your own good. Since the coding is right here, all a lecturer has to do is copy-paste a tiny section into a search engine and they'll find it. Or they could just type in "coding for pi" as well. Believe me, lecturers can have some badass Google-fu skillz.
The main program.
- Minor runtime error in which original main() was declared as 'void' instead of 'int'. While it still works, C doesn't like that.
- Program can now ask for how many decimal places you actually want it to print.
#include "Project Pi(header).h"
long LogBase=4; /* Log10(base) */
/*
** Computation of the constant Pi with arctan relations
*/
int main () {
clock_t endclock, startclock;
long NumDigits, NumArctan;
printf("Please enter number of digits to print pi to: ");
scanf("%ld", &NumDigits);
printf("\n");
long p[10], m[10];
long size=1+NumDigits/LogBase, i;
long *Pi = (long *)malloc(size*sizeof(long));
long *arctan = (long *)malloc(size*sizeof(long));
long *buffer1 = (long *)malloc(size*sizeof(long));
long *buffer2 = (long *)malloc(size*sizeof(long));
startclock = clock();
/*
** Formula used:
**
** Pi/4 = 12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss)
*/
NumArctan = 3;
m[0] = 12; m[1] = 8; m[2] = -5;
p[0] = 18; p[1] = 57; p[2] = 239;
SetToInteger (size, Pi, 0);
/*
** Computation of Pi/4 = Sum(i) [m[i]*arctan(1/p[i])]
*/
for (i=0; i
Mul (size, arctan, abs(m[i]));
if (m[i]>0) Add (size, Pi, arctan);
else Sub (size, Pi, arctan);
}
Mul (size, Pi, 4);
endclock = clock ();
Print (size, Pi); /* Print out of Pi */
printf ("Computation time is : %9.2f seconds\n",
(float)(endclock-startclock)/(float)CLOCKS_PER_SEC );
free (Pi);
free (arctan);
free (buffer1);
free (buffer2);
system ("PAUSE");
}
The functions file.
- Changed the layout and way in which the functions are called/defined. Extra curly brackets basically.
#include "Project Pi(header).h"
const long Base=10000; /* Working base */
long MaxDiv=450; /* about sqrt(2^31/B) */
/*
** Set the big real x to the small integer Integer
*/
void SetToInteger (long n, long *x, long Integer){
long i;
for (i=1; i
x[0] = Integer;
}
}
/*
** Is the big real x equal to zero ?
*/
long IsZero (long n, long *x){
long i;
for (i=0; i
return 0;
return 1;
}
}
}
/*
** Addition of big reals : x += y
** Like school addition with carry management
*/
void Add (long n, long *x, long *y){
long carry=0, i;
for (i=n-1; i>=0; i--){
x[i] += y[i]+carry;
if (x[i]
}
else{
carry = 1;
x[i] -= Base;
}
}
}
/*
** Substraction of big reals : x -= y
** Like school substraction with carry management
** x must be greater than y
*/
void Sub (long n, long *x, long *y){
long i;
for (i=n-1; i>=0; i--){
x[i] -= y[i];
if (x[i]<0){
if (i){
x[i] += Base;
x[i-1]--;
}
}
}
}
/*
** Multiplication of the big real x by the integer q
** x = x*q.
** Like school multiplication with carry management
*/
void Mul (long n, long *x, long q){
long carry=0, xi, i;
for (i=n-1; i>=0; i--){
xi = x[i]*q;
xi += carry;
if (xi>=Base){
carry = xi/Base;
xi -= (carry*Base);
}
else
carry = 0;
x[i] = xi;
}
}
/*
** Division of the big real x by the integer d
** The result is y=x/d.
** Like school division with carry management
** d is limited to MaxDiv*MaxDiv.
*/
void Div (long n, long *x, long d, long *y){
long carry=0, xi, q, i;
for (i=0; i
q = xi/d;
carry = xi-q*d;
y[i] = q;
}
}
/*
** Find the arc cotangent of the integer p = arctan (1/p)
** Result in the big real x (size n)
** buf1 and buf2 are two buffers of size n
*/
void arccot (long p, long n, long *x, long *buf1, long *buf2) {
long p2=p*p, k=3, sign=0;
long *uk=buf1, *vk=buf2;
SetToInteger (n, x, 0);
SetToInteger (n, uk, 1); /* uk = 1/p */
Div (n, uk, p, uk);
Add (n, x, uk); /* x = uk */
while (!IsZero(n, uk)) {
if (p
else {
Div (n, uk, p, uk); /* Two steps for large p (see division) */
Div (n, uk, p, uk);
}
/* uk = u(k-1)/(p^2) */
Div (n, uk, k, vk); /* vk = uk/k */
if (sign) Add (n, x, vk); /* x = x+vk */
else Sub (n, x, vk); /* x = x-vk */
k+=2;
sign = 1-sign;
}
}
/*
** Print the big real x
*/
void Print (long n, long *x) {
long i;
printf ("%d.", x[0]);
for (i=1; i
if (i%25==0) printf ("%8d\n", i*4);
}
printf ("\n");
}
The header file.
- Added function protocols so that the program would work in seperate parts.
- The header file does little, except to make the main and functions files less clunky, and also to eliminate repeating of code.
#include
#include
#include
#include
void SetToInteger (long n, long *x, long Integer);
long IsZero (long n, long *x);
void Add (long n, long *x, long *y);
void Sub (long n, long *x, long *y);
void Mul (long n, long *x, long q);
void Div (long n, long *x, long d, long *y);
void arccot (long p, long n, long *x, long *buf1, long *buf2);
void Print (long n, long *x);
And this following flag counter is simply for my own amusement, to see how many of you guys are actually reading or making use of this stuff. (^_^) Hope it helps.
The original source code was gleaned from the following link (piclassic.c), and was written by Pascal Sebah.
http://numbers.computation.free.fr/Constants/constants.html
This algorithm utilizes the extremely efficient and rapidly converging Machin-like method and formulae to find the value of Pi(π) accurate to up to 10000 decimal places.
When running the program in newer versions of C, however, I found that despite the program working correctly, it printed pi to 10000 decimal places arbitrarily and also seemed to show some runtime errors. I have fixed these, rearranged the program from a single source file to two plus a header to make it neater, changed the function layouts (did old versions of C call them differently?), and also altered it to allow the user to choose the number of digits they would like pi to be printed to.
Original base code however, goes to Pascal Sebah.
If you are using this for an assignment in school/university/etc. then credit both me and Pascal Sebah and for goodness sake change some of the coding or the layout of it so that it's at least some of your own work. At the very least, alter the names of the project, main, functions and header source files. And oh, if it doesn't work, it's because you didn't alter the name of the header file, nor the coding to reflect that.
Remember: Crediting is your friend. I'm not trying to be nasty here, it really is for your own good. Since the coding is right here, all a lecturer has to do is copy-paste a tiny section into a search engine and they'll find it. Or they could just type in "coding for pi" as well. Believe me, lecturers can have some badass Google-fu skillz.
The main program.
- Minor runtime error in which original main() was declared as 'void' instead of 'int'. While it still works, C doesn't like that.
- Program can now ask for how many decimal places you actually want it to print.
#include "Project Pi(header).h"
long LogBase=4; /* Log10(base) */
/*
** Computation of the constant Pi with arctan relations
*/
int main () {
clock_t endclock, startclock;
long NumDigits, NumArctan;
printf("Please enter number of digits to print pi to: ");
scanf("%ld", &NumDigits);
printf("\n");
long p[10], m[10];
long size=1+NumDigits/LogBase, i;
long *Pi = (long *)malloc(size*sizeof(long));
long *arctan = (long *)malloc(size*sizeof(long));
long *buffer1 = (long *)malloc(size*sizeof(long));
long *buffer2 = (long *)malloc(size*sizeof(long));
startclock = clock();
/*
** Formula used:
**
** Pi/4 = 12*arctan(1/18)+8*arctan(1/57)-5*arctan(1/239) (Gauss)
*/
NumArctan = 3;
m[0] = 12; m[1] = 8; m[2] = -5;
p[0] = 18; p[1] = 57; p[2] = 239;
SetToInteger (size, Pi, 0);
/*
** Computation of Pi/4 = Sum(i) [m[i]*arctan(1/p[i])]
*/
for (i=0; i
Mul (size, arctan, abs(m[i]));
if (m[i]>0) Add (size, Pi, arctan);
else Sub (size, Pi, arctan);
}
Mul (size, Pi, 4);
endclock = clock ();
Print (size, Pi); /* Print out of Pi */
printf ("Computation time is : %9.2f seconds\n",
(float)(endclock-startclock)/(float)CLOCKS_PER_SEC );
free (Pi);
free (arctan);
free (buffer1);
free (buffer2);
system ("PAUSE");
}
The functions file.
- Changed the layout and way in which the functions are called/defined. Extra curly brackets basically.
#include "Project Pi(header).h"
const long Base=10000; /* Working base */
long MaxDiv=450; /* about sqrt(2^31/B) */
/*
** Set the big real x to the small integer Integer
*/
void SetToInteger (long n, long *x, long Integer){
long i;
for (i=1; i
x[0] = Integer;
}
}
/*
** Is the big real x equal to zero ?
*/
long IsZero (long n, long *x){
long i;
for (i=0; i
return 0;
return 1;
}
}
}
/*
** Addition of big reals : x += y
** Like school addition with carry management
*/
void Add (long n, long *x, long *y){
long carry=0, i;
for (i=n-1; i>=0; i--){
x[i] += y[i]+carry;
if (x[i]
}
else{
carry = 1;
x[i] -= Base;
}
}
}
/*
** Substraction of big reals : x -= y
** Like school substraction with carry management
** x must be greater than y
*/
void Sub (long n, long *x, long *y){
long i;
for (i=n-1; i>=0; i--){
x[i] -= y[i];
if (x[i]<0){
if (i){
x[i] += Base;
x[i-1]--;
}
}
}
}
/*
** Multiplication of the big real x by the integer q
** x = x*q.
** Like school multiplication with carry management
*/
void Mul (long n, long *x, long q){
long carry=0, xi, i;
for (i=n-1; i>=0; i--){
xi = x[i]*q;
xi += carry;
if (xi>=Base){
carry = xi/Base;
xi -= (carry*Base);
}
else
carry = 0;
x[i] = xi;
}
}
/*
** Division of the big real x by the integer d
** The result is y=x/d.
** Like school division with carry management
** d is limited to MaxDiv*MaxDiv.
*/
void Div (long n, long *x, long d, long *y){
long carry=0, xi, q, i;
for (i=0; i
q = xi/d;
carry = xi-q*d;
y[i] = q;
}
}
/*
** Find the arc cotangent of the integer p = arctan (1/p)
** Result in the big real x (size n)
** buf1 and buf2 are two buffers of size n
*/
void arccot (long p, long n, long *x, long *buf1, long *buf2) {
long p2=p*p, k=3, sign=0;
long *uk=buf1, *vk=buf2;
SetToInteger (n, x, 0);
SetToInteger (n, uk, 1); /* uk = 1/p */
Div (n, uk, p, uk);
Add (n, x, uk); /* x = uk */
while (!IsZero(n, uk)) {
if (p
else {
Div (n, uk, p, uk); /* Two steps for large p (see division) */
Div (n, uk, p, uk);
}
/* uk = u(k-1)/(p^2) */
Div (n, uk, k, vk); /* vk = uk/k */
if (sign) Add (n, x, vk); /* x = x+vk */
else Sub (n, x, vk); /* x = x-vk */
k+=2;
sign = 1-sign;
}
}
/*
** Print the big real x
*/
void Print (long n, long *x) {
long i;
printf ("%d.", x[0]);
for (i=1; i
if (i%25==0) printf ("%8d\n", i*4);
}
printf ("\n");
}
The header file.
- Added function protocols so that the program would work in seperate parts.
- The header file does little, except to make the main and functions files less clunky, and also to eliminate repeating of code.
#include
#include
#include
#include
void SetToInteger (long n, long *x, long Integer);
long IsZero (long n, long *x);
void Add (long n, long *x, long *y);
void Sub (long n, long *x, long *y);
void Mul (long n, long *x, long q);
void Div (long n, long *x, long d, long *y);
void arccot (long p, long n, long *x, long *buf1, long *buf2);
void Print (long n, long *x);
And this following flag counter is simply for my own amusement, to see how many of you guys are actually reading or making use of this stuff. (^_^) Hope it helps.