* add pi calculation

2 files changed, 298 insertions, 1 deletions

diff --git a/src/GNUmakefile.am b/src/GNUmakefile.am
index 83a1a7d..3bb9360 100644
--- a/src/GNUmakefile.am
+++ b/src/GNUmakefile.am
@@ -1,6 +1,8 @@
 bin_PROGRAMS = \
 	float_bench \
-	ff_bench
+	ff_bench \
+	pi_int \
+	pi_fpu
 
 AM_CPPFLAGS = \
 	-I$(top_builddir)
@@ -13,5 +15,12 @@ ff_bench_SOURCES = \
 	ff_bench.c
 ff_bench_LDADD = -lm
 
+pi_int_SOURCES = \
+	pi8.c
+
+pi_fpu_SOURCES = \
+	pi8.c
+pi_fpu_CPPFLAGS = -DUSEFPU
+
 MAINTAINERCLEANFILES = \
 	GNUmakefile.in
diff --git a/src/pi8.c b/src/pi8.c
new file mode 100644
index 0000000..2871a70
--- /dev/null
+++ b/src/pi8.c
@@ -0,0 +1,288 @@
+/* +++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, define:
+
+#define USEFPU     1
+
+  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.
+*/
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+#define SHOWTIME
+
+#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;
+
+
+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 */
+      }
+}
+
+int main(int argc, char *argv[])
+{
+      INDEXER x;
+      time_t T1, T2;
+
+      if (argc != 2)
+      {
+            printf("I need to know how many digits to compute.\n");
+            exit(EXIT_FAILURE);
+      }
+
+      size = atol(argv[1]);
+      if (size <= 0L)
+      {
+            printf("Invalid argument.\n");
+            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))
+      {
+            printf("Unable to allocate enough memory.\n");
+            exit(EXIT_FAILURE);
+      }
+
+      for (x = 0; x < size; x++)
+            pi[x] = 0;
+
+      T1 = time(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
+
+      T2 = time(NULL);
+
+      Print(pi);
+
+#ifdef SHOWTIME
+      printf("\nCalculation time %0.0lf\n", difftime(T2, T1));
+#endif
+
+      return EXIT_SUCCESS;
+}