/* PI calculator by Landon Rabern 1999*/
#include <math.h>
#include <time.h>

int main()
{
  long i,j,q,rmd,rmd2,rmd3,rmd4,max,nl,*pterm,*total;
  float a,b;
  clock_t s,e;
  printf("How many digits: ");
  scanf("%ld",&max);
  pterm=(long *)malloc(max+4);
  total=(long *)malloc(max+4);
  max=max>>2;
  s=clock();
  b=log(5)/log(10000); 
  a=log(0.3125)/log(10000)+2-3*b;
  nl=2;
  pterm[1]=total[1]=3;
  pterm[2]=2000;
  total[2]=1833;

  j=3;
  rmd=870000;
  while (j<=max)
  {
    total[j]=-rmd / 239;
    rmd=rmd % 239;
    j++;
    rmd*=10000;
  }

  i=3;
  do
  {
    j=a+b*i;
    rmd=pterm[j];
    rmd2=0;
    while (j<=max)
    {
      pterm[j]=rmd / 25;
      rmd=10000*(rmd-pterm[j]*25);
      rmd2=rmd2*10000+pterm[j];
      total[j]-=rmd2/i;
      rmd2=rmd2 % i;
      j++;
      if (j<=nl)
        rmd+=pterm[j];
      else if (rmd==0)
      {
        rmd=j;
	while ((rmd2>0) && (j<=max))
	{
	  rmd2=rmd2*10000;
          total[j]-=rmd2 / i;
          rmd2=rmd2 % i;
	  j++;
	}
        j=rmd;
	break;
      }
    }
    nl=--j;
    i+=2;
    j=a+b*i;
    rmd3=pterm[j];
    rmd4=0;
    while (j<=max)
    {
      pterm[j]=rmd3 / 25;
      rmd3=10000*(rmd3-pterm[j]*25);
      rmd4=rmd4*10000+pterm[j];
      total[j]+=rmd4 / i;
      rmd4=rmd4 % i;
      j++;
      if (j<=nl)
        rmd3+=pterm[j];
      else if (rmd3==0)
      {
        rmd3=j;
        while ((rmd4>0) && (j<=max))
	{
          rmd4=rmd4*10000;
          total[j]+=rmd4 / i;
          rmd4=rmd4 % i;
	  j++;
	}
        j=rmd3;
	break;
      }
    }
    nl=--j;
    i+=2;
  } while (i<=max*5.8);
  j=1;
  rmd=4;
  while (j<=max)
  {
    pterm[j]=rmd / 239;
    rmd=rmd % 239;
    j++;
    rmd=rmd*10000;
  }
  b=log(239)/log(10000);
  a=log(59.75)/log(10000)+2-3*b;
  i=3;
  do
  {
    q=i+2;
    j=a+b*i;
    rmd=pterm[j++];
    rmd2=rmd3=rmd4=0;
    while (j<=max)
    {
      rmd=rmd*10000+pterm[j];
      pterm[j]=rmd / 57121;
      rmd-=pterm[j]*57121;
      rmd2=rmd2*10000+pterm[j];
      total[j]+=rmd2/i;
      rmd2=rmd2 % i;

      rmd3=rmd3*10000+pterm[j];
      pterm[j]=rmd3 / 57121;
      rmd3-=pterm[j]*57121;
      rmd4=rmd4*10000+pterm[j];
      total[j]-=rmd4/q;
      rmd4=rmd4 % q;
      j++;
    }
    i+=4;
  } while (i<=max*1.7);
  for (i=max;i>=1;i--)
  {
    if (total[i]<0)
    {
      total[i-1]+=total[i] / 10000-1;
      total[i]=10000+total[i] % 10000;
    }
    if (total[i]>=10000)
    {
      total[i-1]+=total[i] / 10000;
      total[i]=total[i] % 10000;
    }
  }
  e = clock(); 
  printf("%d",3);
  for (j=2;j<=max-1;j++)
  {
    printf("%ld%ld%ld%ld",total[j] / 1000,(total[j] % 1000) / 100,
    ((total[j] % 1000) % 100) / 10,
    ((total[j] % 1000) % 100) % 10);
    printf(" ");
  }
  printf("\n Time = %2.3f seconds\n",(float)(e-s)/CLOCKS_PER_SEC);
  free(total);
  free(pterm);
  return(0);
}