Drexel dragonThe Math ForumDonate to the Math Forum

Ask Dr. Math - Questions and Answers from our Archives
_____________________________________________
Associated Topics || Dr. Math Home || Search Dr. Math
_____________________________________________

A Program to Compute Pi


Date: 06/09/99 at 15:27:58
From: Jan Rembold
Subject: Pi

Hello Dr. Math,

Can you tell me a good algorithm to calculate pi? I found some but 
they didn't work. I want to calculate pi as exactly as possible with 
more than 100 or more than 1000 figures.

Please help me.
Thanks,
Jan


Date: 06/09/99 at 18:01:57
From: Doctor Wilkinson
Subject: Re: Pi

Hi, Jan! 

If you follow the various links from the Pi FAQ at

  http://mathforum.org/dr.math/faq/faq.pi.html   

you can see the first 50,000 digits of pi and you can get some idea of 
the various algorithms that have been used to calculate pi to high 
degrees of accuracy. 

If you want to compute for yourself a lot of digits of pi, you need 
ways to do arithmetic on numbers with a lot of digits, which is a 
problem since the built-in arithmetic in most computers has only 
limited accuracy. Here is a C program that will allow you to compute 
as many digits as you want within reason. If you tried to compute 
100,000,000 decimal places the algorithm (Machin's Formula) is not 
good enough.

Here's the program:

/****************************************/
/* Compute pi to arbitrary precision    */
/* Author Roy Williams February 1994    */
/* Uses Machin's formula...             */
/* pi/4 = 4 arctan(1/5) - arctan(1/239) */
/****************************************/
/* compile with cc -O -o pi pi.c        */
/* run as "pi 1000"                     */
/****************************************/
/* The last few digits may be wrong.......... */

#include <stdio.h>

#define BASE    10000
int nblock;
int *tot;
int *t1;
int *t2;
int *t3;

main(argc, argv)
int argc;
char **argv;
{
        int ndigit = 60;
        if(argc == 2)
                ndigit = atoi(argv[1]);
        else {
                fprintf(stderr, "Usage: %s ndigit\n", argv[0]);
                fprintf(stderr, "(Assuming 60 digits)\n");
        }
        if(ndigit < 20) ndigit = 20;
        nblock = ndigit/4;
        tot = (int *)malloc(nblock*sizeof(int));
        t1 = (int *)malloc(nblock*sizeof(int));
        t2 = (int *)malloc(nblock*sizeof(int));
        t3 = (int *)malloc(nblock*sizeof(int));
        if(!tot || !t1 || !t2 || !t3){
                fprintf(stderr, "Not enough memory\n");
                exit(1);
        }

        arctan(tot, t1, t2, 5, 1);
        mult(tot, 4);
        arctan(t3, t1, t2, 239, 2);
        sub(tot, t3);
        mult(tot, 4);
        print(tot);
}

arctan(result, w1, w2, denom, onestep)
int *result, *w1, *w2, denom, onestep;
{
        int denom2 = denom*denom;
        int k = 1;

        set(result, 1);
        div(result, denom);
        copy(w1, result);

        do{
                if(onestep)
                        div(w1, denom2);
                else {
                        div(w1, denom);
                        div(w1, denom);
                }
                copy(w2, w1);
                div(w2, 2*k+1);
                if(k%2)
                        sub(result, w2);
                else
                        add(result, w2);
                k++;
        } while(!zero(w2));
}

copy(result, from)
int *result, *from;
{
        int i;
        for(i=0; i<nblock; i++)
                result[i] = from[i];
}

zero(result)
int *result;
{
        int i;
        for(i=0; i<nblock; i++)
                if(result[i])
                        return 0;
        return 1;
}

add(result, increm)
int *result, *increm;
{
        int i;
        for(i=nblock-1; i>=0; i--){
                result[i] += increm[i];
                if(result[i] >= BASE){
                        result[i] -= BASE;
                        result[i-1]++;
                }
        }
}

sub(result, decrem)
int *result, *decrem;
{
        int i;
        for(i=nblock-1; i>=0; i--){
                result[i] -= decrem[i];

                if(result[i] < 0){
                        result[i] += BASE;
                        result[i-1]--;
                }
        }
}

mult(result, factor)
int *result, factor;
{
        int i, carry = 0;
        for(i=nblock-1; i>=0; i--){
                result[i] *= factor;
                result[i] += carry;
                carry = result[i]/BASE;
                result[i] %= BASE;
        }
}

div(result, denom)
int *result, denom;
{
        int i, carry = 0;

        for(i=0; i<nblock; i++){
                result[i] += carry*BASE;
                carry = result[i] % denom;
                result[i] /= denom;
        }
}

set(result, rhs)
int *result, rhs;
{
        int i;
        for(i=0; i<nblock; i++)
                result[i] = 0;
        result[0] = rhs;
}

print(result)
int *result;
{
        int i, k;
        char s[10];
        printf("%1d.\n", result[0]);
        for(i=1; i<nblock; i++){
                sprintf(s, "%4d ", result[i]);
                for(k=0; k<5; k++)
                        if(s[k] == ' ') s[k] = '0';
                printf("%c%c%c%c", s[0], s[1], s[2], s[3]);
                if(i%15 == 0) printf("\n");
        }
        printf("\n");
}

This is in classic "K and R" C. 

- Doctor Wilkinson, The Math Forum
  http://mathforum.org/dr.math/   
    
Associated Topics:
High School Calculators, Computers
High School Projects
Middle School Pi

Search the Dr. Math Library:


Find items containing (put spaces between keywords):
 
Click only once for faster results:

[ Choose "whole words" when searching for a word like age.]

all keywords, in any order at least one, that exact phrase
parts of words whole words

Submit your own question to Dr. Math

[Privacy Policy] [Terms of Use]

_____________________________________
Math Forum Home || Math Library || Quick Reference || Math Forum Search
_____________________________________

Ask Dr. MathTM
© 1994-2013 The Math Forum
http://mathforum.org/dr.math/