

Berichten: 1646
Berichten: 98
Locatie: United States
Geslacht: Man
|
Taken from http://nest.cs.uiowa.edu/22C178f00/weblog/hshen/mypi3.c (linked from http://nest.cs.uiowa.edu/22C178f00/weblog/hshen/index.html ). A C program to compute PI to any desired precision. So I win -- I give you infinity digits of PI . (A few ' ' were inserted to prevent emotes from being displayed.)
Program begins below this line.
-------------------------------------------
#include <stdio.h>
#include <gmp.h>
#include <sys/times.h>
#include <sys/time.h>
#include <fcntl.h>
/* to compile: gcc -O2 -o mypi3 mypi3.c -lgmp */
/* to run:
% mypi3 1000
This will give about 1200 digits of PI
*/
unsigned long int ndigit;
#define k1 545140134
#define k2 13591409
#define k3 640320
#define k4 100100025
#define k5 327843840
#define k6 53360
MP_INT s,s1,s2,s3,s4,s5,mypi, temp, temp1, base, an, bn, cn, dn, sa, sb;
main (int argc, char ** argv)
{
int i,j;
struct timeval mytime;
double time1, time2;
struct timezone tzp;
unsigned long int s2_addon;
FILE * fout;
char * numberstr,* kar;
int vijf,tien;
if (argc > 1)
{
ndigit = (unsigned long int) atoi (argv[1]);
}
else
{
ndigit = 40;
}
ndigit = ndigit;
mpz_init (&s);
mpz_init (&s1);
mpz_init (&s2);
mpz_init (&s3);
mpz_init (&s4);
mpz_init (&s5);
mpz_init (&mypi);
mpz_init (&temp);
mpz_init (&temp1);
mpz_init (&base);
mpz_init (&an);
mpz_init (&bn);
mpz_init (&cn);
mpz_init (&dn);
mpz_init (&sa);
mpz_init (&sb);
gettimeofday(&mytime,&tzp);
time1 = mytime.tv_sec + mytime.tv_usec/1000000.0;
/* s1=k6*sqrt(k3) */
mpz_ui_pow_ui(&base,10,ndigit);
mpz_mul_ui(&temp,&base,k3);
mpz_mul(&s1,&temp,&base);
mpz_sqrt(&temp,&s1);
mpz_mul_ui(&s1,&temp,k6);
mpz_set(&s2,&base);
mpz_set_ui(&s3,k2);
mpz_set_ui(&s4,1);
mpz_set_ui(&temp,k4);
mpz_mul_ui(&s5, &temp, k5);
mpz_mul(&s,&s3,&base);
/*
for (i=1; i<ndigit/12; i++)
{
if ((6*i-1)<1600){
s2_addon=(6*i-1)*(6*i-3)*(6*i-5);
mpz_mul_ui(&temp,&s2,s2_addon);
mpz_div_ui(&temp1,&temp,i*i*i);
mpz_div(&s2,&temp1,&s5);
}
else{
mpz_mul_ui(&temp,&s2,6*i-1);
mpz_mul_ui(&s2,&temp,6*i-3);
mpz_mul_ui(&temp,&s2,6*i-5);
mpz_ui_pow_ui(&s2,i,3);
mpz_mul(&temp1,&s2,&s5);
mpz_div(&s2,&temp,&temp1);
}
mpz_add_ui(&temp, &s3, k1);
mpz_set(&s3, &temp);
mpz_mul(&temp1, &s2, &s3);
if((i%2)==1){
mpz_sub(&temp, &s, &temp1);
mpz_set(&s, &temp);
}else{
mpz_add(&temp, &s, &temp1);
mpz_set(&s, &temp);
}
}
*/
for (i=1; i<=ndigit/12/50; i++)
{
mpz_set_ui(&an,0);
mpz_set_ui(&bn,0);
mpz_set_ui(&dn,0);
mpz_set_ui(&sa,1);
mpz_set_ui(&sb,1);
for (j=i*50; j>(i-1)*50; j--)
{
/*cn=(k2+j*k1);*/
mpz_set_ui(&cn, k1);
mpz_mul_ui(&temp, &cn, j);
mpz_add_ui(&cn,&temp,k2);
/* bn=j^3*k4*k5; */
mpz_ui_pow_ui(&temp, j, 3);
mpz_mul(&bn, &temp, &s5);
/* dn=cn*sb-an*dn; */
mpz_mul(&temp,&cn,&sb);
mpz_mul(&temp1, &an, &dn);
mpz_sub(&dn, &temp, &temp1);
/* sb=sb*bn; */
mpz_mul(&temp, &sb, &bn);
mpz_set(&sb, &temp);
/* an=(6*j-1)*(6*j-3)*(6*j-5);*/
mpz_set_ui(&an,6*j-1);
mpz_mul_ui(&temp,&an,6*j-3);
mpz_mul_ui(&an,&temp,6*j-5);
/* sa=sa*an; */
mpz_mul(&temp, &sa, &an);
mpz_set(&sa, &temp);
}
/* dn=an*dn */
mpz_mul(&temp, &an, &dn);
mpz_mul(&dn, &temp, &s2);
/* s=s-dn/sb; */
mpz_div(&temp, &dn, &sb);
mpz_sub(&temp1, &s, &temp);
mpz_set(&s, &temp1);
/* s2=s2*sa/sb; */
mpz_mul(&temp, &s2, &sa);
mpz_div(&s2, &temp, &sb);
}
mpz_mul(&temp,&s1,&base);
mpz_div(&mypi,&temp,&s);
gettimeofday(&mytime,&tzp);
time2 = mytime.tv_sec + mytime.tv_usec/1000000.0;
printf("Time to compute PI: %lf\n\n", time2-time1);
numberstr = NULL;
numberstr = mpz_get_str (numberstr, 10, &mypi);
printf ("%c.\n", numberstr[0]);
vijf = 0;
tien = 0;
kar = numberstr + 1;
while (kar[0] != '\0' ){
putchar (kar[0]);
tien++;
if (tien == 10){
putchar (' ' );
tien = 0;
vijf++;
if (vijf == 5){
putchar ('\n' );
vijf = 0;
}
}
kar++;
}
putchar ('\n' );
free (numberstr);
printf("\n" );
mpz_clear (&s);
mpz_clear (&s1);
mpz_clear (&s2);
mpz_clear (&s3);
mpz_clear (&s4);
mpz_clear (&s5);
mpz_clear (&mypi);
mpz_clear (&temp);
mpz_clear (&temp1);
mpz_clear (&base);
mpz_clear (&an);
mpz_clear (&bn);
mpz_clear (&cn);
mpz_clear (&dn);
mpz_clear (&sa);
mpz_clear (&sb);
return;
}
__________________________
SkyLords Head Programmer
Spelled: I I R I (not irii, irri, or iri).
Force of nature.
|