Source Code in C for pi

[Geoffrey Faivre-Malloy]

I'd like to get source code in C for pi (hopefully this will run on a
unix box). Any pointers? Also, any pointers to prime number source
code would be appreciated. Please E-Mail me as I don't frequent this
newsgroup. Thanks in advance,

Geoffrey

Geo...@crl.com

-------------------------------Sig? What Sig?-------------------------------

[Aephraim M. Steinberg]

In article <2lt9tm$5...@crl2.crl.com>,

Geoffrey Faivre-Malloy <geof...@crl.com> wrote:
>I'd like to get source code in C for pi (hopefully this will run on a
>unix box). Any pointers? Also, any pointers to prime number source

Here's a two-line program someone in my group ran across, which I have
never succeeded in understanding. It's fast. Also not too hard to
modify it to get more digits. (By the same token, not too hard to
modify it to give gibberish, so long as you don't understand the
algorithm! I'm convinced that it already gives pi to a much better
approximation than the degree to which our universe is Euclidean
anyway :-)

int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;for(;d=0,g=c*2;c
-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}

A slightly de-obfuscated version follows, but we don't understand
it either...

#define SCALE 10000
#define MAXARR 2800
#define ARRINIT 2000

main()
{
int i, j;
int carry = 0;
int arr[MAXARR+1];

for (i = 0; i <= MAXARR; ++i)
arr[i] = ARRINIT;
for (i = MAXARR; i; i -= 14) {
int sum = 0;
for (j = i; j > 0; --j) {
sum = sum*j + SCALE*arr[j];
arr[j] = sum % (j*2-1);
sum /= (j*2-1);
}
printf("%04d", carry + sum/SCALE);
carry = sum % SCALE;
}
}

--
Aephraim M. Steinberg | "It's no exaggeration to say
UCB Physics | the undecideds could go one way
aeph...@physics.berkeley.edu | or the other."
| -- George Bush

[David E. Brown]

aephraim@physics5 (Aephraim M. Steinberg) writes:

>In article <2lt9tm$5...@crl2.crl.com>,
>Geoffrey Faivre-Malloy <geof...@crl.com> wrote:
>>I'd like to get source code in C for pi (hopefully this will run on a
>>unix box). Any pointers? Also, any pointers to prime number source

>Here's a two-line program someone in my group ran across, which I have

[elided]

How about 4*atan(1.0)?

[Dik T. Winter]

In article <2m2g0t$g...@agate.berkeley.edu> aephraim@physics5 (Aephraim M. Steinberg) writes:
> Here's a two-line program someone in my group ran across, which I have
> never succeeded in understanding.

It is not so difficult I think, but then, I wrote it. It appears also to
give the best compression for the first 800 digits of pi ;-).

> Also not too hard to
> modify it to get more digits.

It is not hard to modify it to get more digits. It is hard to modify it
to get more correct digits. The integers 10000, 2800, 2801, 14 and 4
are strongly connected.
--
dik t. winter, cwi, kruislaan 413, 1098 sj amsterdam, nederland, +31205924098
home: bovenover 215, 1025 jn amsterdam, nederland; e-mail: d...@cwi.nl

[Kerry Cianos]

Geoffrey,

Here is some C code to compute Pi for a specified number of
digits, assuming you have enough memory :) I found this in one
of the feeds a year or two ago. Enjoy!

Kerry Cianos
Senior Software Engineer, Purple Duck Software
eur...@crl.com

--- cut ---

#include <stdio.h>

long kf,ks,*mf,*ms;
long cnt,n,temp,nd;
long i;
long col,col1;
long loc,stor[21];
FILE *fp;

shift(l1,l2,lp,lmod)
long *l1,*l2,lp,lmod;
{
long k;
k=((*l2)>0 ? (*l2)/lmod : -(-(*l2)/lmod)-1);
*l2 -= k*lmod;
*l1 += k*lp;
}

yprint(m)
long m;
{

if (cnt < n)
{
if (++col==11)
{
col = 1;
if (++col1==6)
{
col1=0;
fprintf(fp,"\n");
printf("\n");
fprintf(fp,"%4ld",m%10);
printf("%4ld",m%10);
}
else {
fprintf(fp,"%3ld",m%10);
printf("%3ld",m%10);
}
}
else {
fprintf(fp,"%ld",m);
printf("%ld",m);
}
cnt++;
}
}

xprint(m)
long m;
{
long ii,wk,wk1;

if (m<8)
{
for (ii=1; ii<=loc;)
yprint(stor[ii++]);
loc=0;
}
else
if (m>9)
{
wk = m/10;
m %= 10;
for (wk1=loc; wk1>=1;wk1--)
{
wk += stor[wk1];
stor[wk1] = wk % 10;
wk /= 10;
}
}
stor[++loc] = m;
}

main(argc,argv)
int argc;
char *argv[];
{
int i=0;

stor[i++] = 0;
if (argc < 2)
{
fprintf(stderr,"Format is:\n\tlongpi <# of places>\n");
exit();
}

n = atoi(argv[1]);

mf = (long *)calloc((int)n+3, sizeof(long));
if (mf==0)
{
fprintf(stderr,"Memory allocation failure [mf]\n");
exit();
}

ms = (long *)calloc((int)n+3, sizeof(long));
if (ms==0)
{
fprintf(stderr,"Memory allocation failure [ms]\n");
exit();
}

printf("Memory allocated.\n");

fp = fopen("pi.out","w");
fprintf(fp,"\nThe following is an approximation of PI to %ld \
digits\n",(long)n);
printf("\nThe following is an approximation of PI to %ld \
digits\n",(long)n);


cnt = 0;
kf = 25;
ks = 57121L;
mf[1] = 1;
for (i=2; i<=n; i+=2)
{
mf[i] = -16;
mf[i+1] = 16;
}

for (i=1; i<=n; i+=2)
{
ms[i] = -4;
ms[i+1] = 4;
}

fprintf(fp,"\n 3.");
printf("\n 3.");

while (cnt < n)
{
for (i=0; ++i<=n-cnt;)
{
mf[i] *= 10;
ms[i] *= 10;
}
for (i=n-cnt+1; --i>=2;)
{
temp = 2*i-1;
shift(&mf[i-1],&mf[i],temp-2,temp*kf);
shift(&ms[i-1],&ms[i],temp-2,temp*ks);
}
nd = 0;
shift(&nd,&mf[1],1L,5L);
shift(&nd,&ms[1],1L,239L);
xprint(nd);
}
printf("\n\nCalculations Completed!\n");
}

[Dik T. Winter]

In article <CMoDs...@cwi.nl> d...@cwi.nl (Dik T. Winter) writes:
> In article <2m2g0t$g...@agate.berkeley.edu> aephraim@physics5 (Aephraim M. Steinberg) writes:
> > Here's a two-line program someone in my group ran across, which I have
> > never succeeded in understanding.
>
> It is not so difficult I think, but then, I wrote it. It appears also to
> give the best compression for the first 800 digits of pi ;-).
>

Within a few hours of posting this I have already received quite a few
requests for explanation. Apart from the two level C obfuscation there
is also an algorithmic obfuscation. I will show how it is done by
example; it is for the reader what is actually happening (and he may
use the unobfuscated source that was posted for that).

Let's try to print 'e' to a large number of decimals. We can successively
calculate '1/n!' for increasing values of 'n' and add them together, but
that is very inefficient. It is better to represent 'e' in a mixed-number
base. The zeroth digit should be multiplied by '1/0!', the first by '1/1!',
the second by '1/2!' and in general the n-th by '1/n!'. To get the next
digit, we multiply by 10, cater for carries and extract the integer part.
In this case, because the zeroth and the first digit both should be
multiplied by 1, we collapse them and start with a first digit of 2.
We get the following C program:
int digits[1000], i, d, carry;

main()
{
for(i = 1; i < 1000; i++) digits[i] = 1; /* the digits in the rep. */
digits[1] = 2; /* collapse zeroth and first digit */
while (1) { /* infinite loop; have a guess when it starts to
produce wrong digits! */
putchar(digits[1] + '0'); digits[1] = 0; /* print digit and clear */
carry = 0; /* initialize carry */
for(i = 999; i > 1; i--) { /* loop down over the representation */
d = digits[i] * 10 + carry;
digits[i] = d % i;
carry = d / i;
}
digits[1] = carry;
}
}

Exercise: Why do the carries work this way, what represention and series
were used for pi?

You may want to print a newline on occasion.