/*
 * Compute many digits of pi with Mach's formula:
 * pi/4 = 4 arctan(1/5) - arctan(1/239)
 * arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + -...
 *
 * Give desired number of digits as argument,
 * add "-n" to output only series lengths.
 *
 */

#include <stdio.h>

#define DIGITSIZE 4	/* # of decimal digits in Digit */
#define BASE 10000 	/*  10^DIGITSIZE */

typedef unsigned short Digit;   	/* must hold 2*BASE+1    */
typedef unsigned long LongDigit;	/* must hold 239^2*BASE  */


void PrintBig(Big, size)
Digit *Big;
unsigned size;
{
  int i, k;
  char s[5], *j;
  printf("\n3.");
  for (i = 0, k=0; i < size; i++) {
    sprintf(s,"%0*u", DIGITSIZE, Big[i]);
    for (j = s; *j; j++) {
	putchar(*j);
	if (++k % 5 == 0) {
	   putchar(' ');
	   if (k % 50 == 0) {
	      printf("\n  ");
	      k = 0;
	   }
	}
    }
  }
  if (k)
     printf("\n");
  printf("\n");
}


void ClearBig(Big, size)
Digit *Big;
unsigned size;
{
  Digit *i;
  for (i = Big; i < &Big[size]; i++) {
    *i = 0;
  }
}


void BigMult(BigFactor,SmallFactor,carry,size)
     Digit *BigFactor;
     Digit SmallFactor;
     Digit carry;
     unsigned size;
{
  Digit *i;
  LongDigit t;
  for (i = &BigFactor[size]; --i >= BigFactor; ) {
    t = (LongDigit)*i * SmallFactor + carry;
    carry = t / BASE;
    *i = t % BASE;
  }
}


Digit *BigDiv(BigDividend,SmallDivisor,rem,size)
     Digit *BigDividend;
     Digit SmallDivisor;
     Digit *rem;
     unsigned size;
{
  Digit *i, *last = 0;
  LongDigit carry = 0, tmp;
  for (i = BigDividend; i < &BigDividend[size]; i++) {
    tmp = *i + carry;
    carry = (tmp % SmallDivisor) * BASE;
    *i = tmp / SmallDivisor;
    if (*i) {
      last = i;
    }
  }
  *rem = carry / BASE;
  return last;
}


void BigAdd(BigSum,BigTerm,size)
     Digit *BigSum;
     Digit *BigTerm;
     unsigned size;
{
  Digit *i, *j, carry = 0;
  for (i = &BigSum[size-1], j = &BigTerm[size-1];
       i >= BigSum;
       i--, j--) {
    *i = *i + *j + carry;
    carry = *i / BASE;
    *i %= BASE;
  }
}


void BigSub(BigSum,BigTerm,size)
     Digit *BigSum;
     Digit *BigTerm;
     unsigned size;
{
  Digit *i, *j, borrow = 0;
  for (i = &BigSum[size-1], j = &BigTerm[size-1];
       i >= BigSum;
       i--, j--) {
    *i += BASE;
    *i = *i - *j - borrow;
    borrow = (*i < BASE);
    *i %= BASE;
  }
}


Digit AddArcTan(BigSum,n,BigTerm,size)
     Digit *BigSum;
     int n;
     Digit *BigTerm;
     unsigned size;
{
  Digit n2 = (Digit)n*n, divisor = 1, rem;

  ClearBig(BigTerm,size);
  BigTerm[0] = BASE;
  if (n < 0) {
      BigDiv(BigTerm, (Digit)-n, &rem, size);
      BigSub(BigSum,BigTerm,size);
      divisor += 2;
    } else {
      BigDiv(BigTerm, (Digit)n, &rem, size);
      BigAdd(BigSum,BigTerm,size);
      divisor += 2;
      BigDiv(BigTerm, n2, &rem, size);
      BigDiv(BigTerm, divisor, &rem, size);
      BigSub(BigSum,BigTerm,size);
      BigMult(BigTerm,divisor,rem,size);
      divisor += 2;
    }

  while (BigDiv(BigTerm, n2, &rem, size)) {
      BigDiv(BigTerm, divisor, &rem, size);
      BigAdd(BigSum,BigTerm,size);
      BigMult(BigTerm,divisor,rem,size);
      divisor += 2;
      BigDiv(BigTerm, n2, &rem, size);
      BigDiv(BigTerm, divisor, &rem, size);
      BigSub(BigSum,BigTerm,size);
      BigMult(BigTerm,divisor,rem,size);
      divisor += 2;
  }
  return divisor;
}


main(argc,argv)
     int argc;
     char **argv;
{
  int decimals, silent;
  unsigned size;
  static Digit *BigSum, *BigTerm;

  if (argc < 2 || 1 != sscanf(argv[1], "%d", &decimals)) {
    fprintf(stderr,"usage:  pi n\n");
    exit(1);
  }
  silent = argc > 2 && !strcmp(argv[2],"-n");
  size = (decimals + DIGITSIZE -1) / DIGITSIZE;
  if (NULL == (BigSum = (Digit *)malloc(size * sizeof(Digit))) ||
      NULL == (BigTerm = (Digit *)malloc(size * sizeof(Digit)))) {
    fprintf(stderr,"insufficient memory\n");
    exit(2);
  }
  ClearBig(BigSum,size);
  printf("last divisor in atan(1/5) series: %u\n",
	 AddArcTan(BigSum,5,BigTerm,size));
  BigMult(BigSum,(Digit)4,(Digit)0,size);
  printf("last divisor in atan(1/239) series: %u\n",
	 AddArcTan(BigSum,-239,BigTerm,size));
  BigMult(BigSum,(Digit)4,(Digit)0,size);
  if (!silent)
    PrintBig(BigSum,size);
  return 0;
}