path: root/src/pi8.c

/* +++Date last modified: 01-Oct-1996 */

/*
  pi8.c  Sept 9, 1996.

  This program is placed into the public domain by its author,
  Carey Bloodworth.

  This pi program can calculate pi with either integers, or doubles.
  It can also use a variety of formulas.

  This code isn't really optimized because it has to work with both the
  long integer version and the double version, and several formulas.
  Compromises had to be made in several places.

  When compiling, you can chose to use the FPU or the integer version.
  By default, it will chose to work only in integers.  If you want to
  use the FPU, compile with additional flag -DUSEFPU

  You have a choice of formulas.  By default, it will use the Machin
  formula of:  pi=16arctan(1/5)-4arctan(1/239)

  You could chose to use one of the other formulas.

  for pi=8arctan(1/3)+4arctan(1/7)
#define ARC3  1
  for pi=12arctan(1/4)+4arctan(1/20)+4arctan(1/1985)
#define ARC4  1
  for pi=16arctan(1/5)-4arctan(1/70)+4arctan(1/99)
#define ARC5  1
  for pi=32arctan(1/10)-4arctan(1/239)-16arctan(1/515)
#define ARC10 1

  Or, of course, you could define it on the compile command line with
  the -D switch.

  Timings were done on a Cyrix 486/66, with the slow Turbo C++ v3.0
          1,000 2,000 3,000 4,000  1,000F 2,000F 3,000F 4,000F
  Machin      4    18    45    86      1      5     11     20
  Arc3        6    29    74   140      2      9     20     35
  Arc4        5    24    64   116      2      7     16     29
  Arc5        6    26    65   123      1      6     15     26
  Arc10       4    19    46    86      1      5     11     19

  All of the combinations above were verified to their indicated
  precision and in each case, only the last few digits were wrong,
  which is a normal round off / truncation error.

  Better compilers will of course result in faster computations,
  but the ratios should be the same.  When I used GCC 2.6.3 for
  DOS, I computed 4,000 digits with with the Machin formula and
  the FPU in 8 seconds.  The integer version took 17 seconds.

  I also used the FPU GCC version to compute 65,536 digits of
  pi and verified them against the Gutenberg PIMIL10.TXT, and
  only the last 4 digits were incorrect.  The computations took
  33 minutes and 54 seconds.
*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#include <unistd.h>

#define SHOWTIME

#define DEFAULT_DIGIT_COUNT 2000

#if defined USEFPU

#define BASE   1000000000L
#define BASEDIGITS 9

typedef long int SHORT;
typedef double   LONG;
#else

#define BASE   100
#define BASEDIGITS 2

typedef unsigned char SHORT;
typedef long int      LONG;
#endif

typedef long int INDEXER;


SHORT *pi, *powers, *term;
INDEXER size = DEFAULT_DIGIT_COUNT;


void OutDig(int dig)
{
      static int printed = 0;

      putchar(dig + '0');
      printed++;
      if ((printed%1000) == 0)
      {
            printed = 0;
            printf("\n\n\n");
      }
      if ((printed%50) == 0)
            printf("\n");
      else if ((printed%10) == 0)
            putchar(' ');
}


void PrintShort(SHORT num)
{
      int x;
      int digits[BASEDIGITS + 1];

      for (x = 0; x < BASEDIGITS; x++)
      {
            digits[x] = num % 10;
            num /= 10;
      }
      for (x = BASEDIGITS - 1; x >= 0; x--)
            OutDig(digits[x]);
}

void Print(SHORT *num)
{
      INDEXER x;

      printf("\nPI = 3.\n");
      for (x = 1; x < size; x++)
            PrintShort(num[x]);
      printf("\n");
}

void arctan(int multiplier, int denom, int sign)
{
      INDEXER x;
      LONG remain, temp, divisor, denom2;
      SHORT NotZero = 1;
      INDEXER adv;

      for (x = 0; x < size; x++)
            powers[x] = 0;

      divisor = 1;
      denom2 = (LONG)denom;denom2 *= denom2;
      adv = 0;

      remain = (LONG)multiplier * denom;
      while (NotZero)
      {
            for (x = adv; x < size; x++)
            {
                  temp = (LONG)powers[x] + remain;
                  powers[x] = (SHORT)(temp / denom2);
                  remain = (temp - (denom2 * (LONG)powers[x])) * BASE;
            }

            remain = 0;
            for (x = adv; x < size; x++)
            {
                  temp = (LONG)powers[x] + remain;
                  term[x] = (SHORT)(temp / divisor);
                  remain = (temp - (divisor * (LONG)term[x])) * BASE;
            }
            remain = 0;

            if (sign > 0)
            {
                  LONG carry, sum;

                  carry = 0;
                  for (x = size - 1; x >=0; x--)
                  {
                        sum = (LONG)pi[x] + (LONG)term[x] + carry;
                        carry = 0;
                        if (sum >= BASE)
                        {
                              carry = 1;
                              sum -= BASE;
                        }
                        pi[x] = (SHORT)sum;
                  }
            }
            else
            {
                  LONG borrow, sum;

                  borrow = 0;
                  for (x = size - 1; x >= 0; x--)
                  {
                        sum = (LONG)pi[x] - (LONG)term[x] - borrow;
                        borrow = 0;
                        if (sum < 0)
                        {
                              borrow = 1;
                              sum += BASE;
                        }
                        pi[x] = (SHORT)sum;
                  }
            }

            sign = -sign;
            divisor += 2;
            NotZero = 0;
            for (x = adv; x < size; x++)
            {
                  if (powers[x])
                  {
                        NotZero = 1;
                        break;
                  }
            }

            if (NotZero)
            {
                  while (powers[adv] == 0)
                        adv++;
            }
            /* We can skip ones that are already 0 */
      }
}

static void usage(const char *prg, FILE *out)
{
	fprintf(out, "usage: %s [options]\n", prg);
	fprintf(out, " options are:\n");
	fprintf(out, " '-c count' digits to compute. Default is %ld, should not be 0.\n", size);
	fprintf(out, " '-q' be quiet, do not emit pi\n");
	fprintf(out, " '-v' print program version and exit\n");
	fprintf(out, " '-h' print this help and exit\n");
}

int main(int argc, char *argv[])
{
	INDEXER x;
	int quiet = 0, c;
	struct timeval start, stop, diff;

	/* handle command line options first */
	while (1) {
		c = getopt(argc, argv, "c:qhv");
		if (c == -1)
			break;

		switch (c) {
		case 'c':
			size = atoi(optarg);	/* digits to compute */
			break;
		case 'q':
			quiet = 1;
			break;
		case 'h':
			usage(argv[0], stdout);
			exit(EXIT_SUCCESS);
		case 'v':
			printf("%s %s\n", PACKAGE_NAME, PACKAGE_VERSION);
			exit(EXIT_SUCCESS);
		}
	}

	if (! size) {
		usage(argv[0], stderr);
		exit(EXIT_FAILURE);
	}

	size = ((size + BASEDIGITS - 1) / BASEDIGITS) + 1;

	pi = malloc(sizeof(SHORT) * size);
	powers = malloc(sizeof(SHORT) * size);
	term = malloc(sizeof(SHORT) * size);
	
	if ((pi == NULL) || (powers == NULL) || (term == NULL)) {
		fprintf(stderr, "Unable to allocate enough memory.\n");
		exit(EXIT_FAILURE);
	}
	
	for (x = 0; x < size; x++)
		pi[x] = 0;
	
	gettimeofday(&start, NULL);

#if defined ARC3
	arctan(8, 3, 1);
	arctan(4, 7, 1);
#elif defined ARC5
	arctan(16, 5, 1);
	arctan(4, 70, -1);
	arctan(4, 99, 1);
#elif defined  ARC4
	arctan(12, 4, 1);
	arctan(4, 20, 1);
	arctan(4, 1985, 1);
#elif defined  ARC10
	arctan(32, 10, 1);
	arctan(4, 239, -1);
	arctan(16, 515, -1);
#else
      /* Machin formula */
	arctan(16, 5, 1);
	arctan(4, 239, -1);
#endif

	gettimeofday(&stop, NULL);
	if (!quiet)
		Print(pi);
	
	timersub(&stop, &start, &diff);
	printf("\nTime (s): %ld.%ld\n", diff.tv_sec, diff.tv_usec);
	
	return EXIT_SUCCESS;
}