calc_pi.c

计算PI的值的示例代码。仅供参考.无版权责任。

#include "stdio.h"
#include "stdlib.h"
#include "unistd.h"

int _profileDemoTID = 0;

/*
int numDigits = 10;
*/

int numDigits = 1000;


void calcPI(int ndigit, int printResult);

/*************************************************************************
 *
 * ProfileDemoMain - main routine for calculating PI
 *
 */
int main()
{
    calcPI(numDigits, 1);
    while(1) {
	calcPI(numDigits, 0);
	usleep(1);
    }	
    return 0;
}

/*************************************************************************
 *

   Code used to calculate PI to an arbitrary precision (memory
   limited).

   Modification History 
   07may02, rti ANSIfied

   Derived from:
   Compute pi to arbitrary precision    
   Author Roy Williams February 1994    
   [e-mail]               
   Uses Machin's formula...            
   pi/4 = 4 arctan(1/5) - arctan(1/239)

   The last few digits may be wrong.......... 
*/


#define BASE    10000

void Arctan(int *result, int denom, int onestep, int nblock);
void Copy(int *result, int *from, int nblock);
int  Zero(int *result, int nblock);
void Sub(int *result, int *decrem, int nblock);
void Mult(int *result, int factor, int nblock);
void Div(int *result, int denom, int nblock);
void Set(int *result, int rhs, int nblock);
void Print(int *result, int nblock);
void Add (int *result, int *increm, int nblock);

void calcPI(int ndigit, int printResult)
{
    int nblock;
    int *tot, *tmp;

    if (ndigit < 20) ndigit = 20;
    if (ndigit > 100000) {
	printf("Not allowed to calculate over 100,000 digits.\n");
	return;
    }

    nblock = ndigit/4;

    tot = (int *)malloc(nblock*sizeof(int));
    tmp = (int *)malloc(nblock*sizeof(int));

    Arctan(tot, 5, 1, nblock);
    Mult(tot, 4, nblock);
    Arctan(tmp, 239, 2, nblock);
    Sub(tot, tmp, nblock);
    Mult(tot, 4, nblock);

    if (printResult) {
       Print(tot, nblock); 
    }

    free(tot);
    free(tmp);
}


void Arctan(int *result, int denom, int onestep, int nblock)
{
    int denom2 = denom*denom;
    int *w1, *w2;
    int k = 1;
    
    w1 = (int *)malloc(nblock*sizeof(int));

    Set(result, 1, nblock);
    Div(result, denom, nblock);
    Copy(w1, result, nblock);
    
    while(1) {
	if(onestep)
	    Div(w1, denom2, nblock);
	else {
	    Div(w1, denom, nblock);
	    Div(w1, denom, nblock);
	}
	w2 = (int *)malloc(nblock*sizeof(int));
	Copy(w2, w1, nblock);
	Div(w2, 2*k+1, nblock);

	if(k%2)
	    Sub(result, w2, nblock);
	else
	    Add(result, w2, nblock);
	k++;

	if (Zero(w2, nblock)) {
	    free(w2);
	    break;
	} else {
	    free(w2);
	}
    };

    free(w1);
}

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

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

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

void Sub(int *result, int *decrem, int nblock)
{
    int i;
    for(i=nblock-1; i>=0; i--){
	result[i] -= decrem[i];
	if(result[i] < 0){
	    result[i] += BASE;
	    result[i-1]--;
	}
    }
}

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

void Div(int *result, int denom, int nblock)
{
    int carry = 0; 
    int i;

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

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

void Print(int *result, int nblock)
{
    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");
}