/*
 * Compute lots of digits on 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 + -...
 *
 * Usable for up to 100000 digits or so before getting too slow
 * (depending on your computer speed of course)
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifndef DIGITSIZE
#define DIGITSIZE 4
#endif
#undef BASE

/* In each case,
   BASE      = 10^DIGITSIZE
   Digit     = smallest type that can hold max(2*BASE+1,239^2)
               (239^2 comes from n2 in AddArcTan and means char will never do)
   LongDigit = smallest type that can hold 239^2*BASE

   Optimal DIGITSIZE values are 4, 9 and 14,
   being the largest that short, long and long long can handle,
   respectively.   
   This version can handle all values from 1 to 14.
 */

#if DIGITSIZE == 1
#define BASE 10
typedef unsigned short Digit;
typedef unsigned long LongDigit;
#define DigitFmt "u"

#elif DIGITSIZE == 2
#define BASE 100
typedef unsigned short Digit;
typedef unsigned long LongDigit;
#define DigitFmt "u"

#elif DIGITSIZE == 3
#define BASE 1000
typedef unsigned short Digit;
typedef unsigned long LongDigit;
#define DigitFmt "u"

#elif DIGITSIZE == 4
#define BASE 10000 	                /* 10^DIGITSIZE */
typedef unsigned short Digit;   	/* must hold 2*BASE+1    */
typedef unsigned long LongDigit;	/* must hold 239^2*BASE  */
#define DigitFmt "u"

#elif DIGITSIZE == 5
#define BASE 100000
typedef unsigned long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "lu"

#elif DIGITSIZE == 6
#define BASE 1000000
typedef unsigned long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "lu"

#elif DIGITSIZE == 7
#define BASE 10000000
typedef unsigned long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "lu"

#elif DIGITSIZE == 8
#define BASE 100000000
typedef unsigned long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "lu"

#elif DIGITSIZE == 9
#define BASE 1000000000
typedef unsigned long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "lu"

#elif DIGITSIZE == 10
#define BASE 10000000000
typedef unsigned long long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "llu"

#elif DIGITSIZE == 11
#define BASE 100000000000
typedef unsigned long long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "llu"

#elif DIGITSIZE == 12
#define BASE 1000000000000
typedef unsigned long long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "llu"

#elif DIGITSIZE == 13
#define BASE 10000000000000
typedef unsigned long long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "llu"

#elif DIGITSIZE == 14
#define BASE 100000000000000
typedef unsigned long long Digit;
typedef unsigned long long LongDigit;
#define DigitFmt "llu"

#endif

#ifndef BASE
#error DIGITSIZE must be between an integer 1...14.
#else

typedef unsigned Index;


void PrintBig(Digit *Big, Index size)
{
  int i, k;
  char s[DIGITSIZE+1], *j;
  printf("\n3.");
  for (i = 0, k=0; i < size; i++) {
    sprintf(s,"%0*"DigitFmt, DIGITSIZE, (Digit)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(Digit *Big, Index size)
{
  Digit *i;
  for (i = Big; i < &Big[size]; i++) {
    *i = 0;
  }
}


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


Digit *BigDiv(Digit *BigDividend, Digit SmallDivisor, Digit *rem, Index 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(Digit *BigSum, Digit *BigTerm, Index 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(Digit *BigSum, Digit *BigTerm, Index 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(Digit *BigSum, int n, Digit *BigTerm, Index 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;
}

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

  if (argc < 2 || 1 != sscanf(argv[1], "%d", &decimals)) {
    fprintf(stderr,"usage:  pi n [-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: %"DigitFmt"\n",
	 AddArcTan(BigSum,5,BigTerm,size));
  BigMult(BigSum,(Digit)4,(Digit)0,size);
  printf("last divisor in atan(1/239) series: %"DigitFmt"\n",
	 AddArcTan(BigSum,-239,BigTerm,size));
  BigMult(BigSum,(Digit)4,(Digit)0,size);
  if (!silent)
    PrintBig(BigSum,size);
  return 0;
}

#endif