home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / snip9707.zip / PI8.C < prev    next >
C/C++ Source or Header  |  1997-07-05  |  8KB  |  289 lines
  1. /* +++Date last modified: 05-Jul-1997 */
  2.  
  3. /*
  4.   pi8.c  Sept 9, 1996.
  5.  
  6.   This program is placed into the public domain by its author,
  7.   Carey Bloodworth.
  8.  
  9.   This pi program can calculate pi with either integers, or doubles.
  10.   It can also use a variety of formulas.
  11.  
  12.   This code isn't really optimized because it has to work with both the
  13.   long integer version and the double version, and several formulas.
  14.   Compromises had to be made in several places.
  15.  
  16.   When compiling, you can chose to use the FPU or the integer version.
  17.   By default, it will chose to work only in integers.  If you want to
  18.   use the FPU, define:
  19.  
  20. #define USEFPU     1
  21.  
  22.   You have a choice of formulas.  By default, it will use the Machin
  23.   formula of:  pi=16arctan(1/5)-4arctan(1/239)
  24.  
  25.   You could chose to use one of the other formulas.
  26.  
  27.   for pi=8arctan(1/3)+4arctan(1/7)
  28. #define ARC3  1
  29.   for pi=12arctan(1/4)+4arctan(1/20)+4arctan(1/1985)
  30. #define ARC4  1
  31.   for pi=16arctan(1/5)-4arctan(1/70)+4arctan(1/99)
  32. #define ARC5  1
  33.   for pi=32arctan(1/10)-4arctan(1/239)-16arctan(1/515)
  34. #define ARC10 1
  35.  
  36.   Or, of course, you could define it on the compile command line with
  37.   the -D switch.
  38.  
  39.   Timings were done on a Cyrix 486/66, with the slow Turbo C++ v3.0
  40.           1,000 2,000 3,000 4,000  1,000F 2,000F 3,000F 4,000F
  41.   Machin      4    18    45    86      1      5     11     20
  42.   Arc3        6    29    74   140      2      9     20     35
  43.   Arc4        5    24    64   116      2      7     16     29
  44.   Arc5        6    26    65   123      1      6     15     26
  45.   Arc10       4    19    46    86      1      5     11     19
  46.  
  47.   All of the combinations above were verified to their indicated
  48.   precision and in each case, only the last few digits were wrong,
  49.   which is a normal round off / truncation error.
  50.  
  51.   Better compilers will of course result in faster computations,
  52.   but the ratios should be the same.  When I used GCC 2.6.3 for
  53.   DOS, I computed 4,000 digits with with the Machin formula and
  54.   the FPU in 8 seconds.  The integer version took 17 seconds.
  55.  
  56.   I also used the FPU GCC version to compute 65,536 digits of
  57.   pi and verified them against the Gutenberg PIMIL10.TXT, and
  58.   only the last 4 digits were incorrect.  The computations took
  59.   33 minutes and 54 seconds.
  60. */
  61.  
  62. #include <stdio.h>
  63. #include <stdlib.h>
  64. #include <time.h>
  65.  
  66. #define SHOWTIME
  67. #define USEFPU
  68.  
  69. #if defined USEFPU
  70.  
  71. #define BASE   1000000000L
  72. #define BASEDIGITS 9
  73.  
  74. typedef long int SHORT;
  75. typedef double   LONG;
  76. #else
  77.  
  78. #define BASE   100
  79. #define BASEDIGITS 2
  80.  
  81. typedef unsigned char SHORT;
  82. typedef long int      LONG;
  83. #endif
  84.  
  85. typedef long int INDEXER;
  86.  
  87.  
  88. SHORT *pi, *powers, *term;
  89. INDEXER size;
  90.  
  91.  
  92. void OutDig(int dig)
  93. {
  94.       static int printed = 0;
  95.  
  96.       putchar(dig + '0');
  97.       printed++;
  98.       if ((printed%1000) == 0)
  99.       {
  100.             printed = 0;
  101.             printf("\n\n\n");
  102.       }
  103.       if ((printed%50) == 0)
  104.             printf("\n");
  105.       else if ((printed%10) == 0)
  106.             putchar(' ');
  107. }
  108.  
  109.  
  110. void PrintShort(SHORT num)
  111. {
  112.       int x;
  113.       int digits[BASEDIGITS + 1];
  114.  
  115.       for (x = 0; x < BASEDIGITS; x++)
  116.       {
  117.             digits[x] = num % 10;
  118.             num /= 10;
  119.       }
  120.       for (x = BASEDIGITS - 1; x >= 0; x--)
  121.             OutDig(digits[x]);
  122. }
  123.  
  124. void Print(SHORT *num)
  125. {
  126.       INDEXER x;
  127.  
  128.       printf("\nPI = 3.\n");
  129.       for (x = 1; x < size; x++)
  130.             PrintShort(num[x]);
  131.       printf("\n");
  132. }
  133.  
  134. void arctan(int multiplier, int denom, int sign)
  135. {
  136.       INDEXER x;
  137.       LONG remain, temp, divisor, denom2;
  138.       SHORT NotZero = 1;
  139.       INDEXER adv;
  140.  
  141.       for (x = 0; x < size; x++)
  142.             powers[x] = 0;
  143.  
  144.       divisor = 1;
  145.       denom2 = (LONG)denom;denom2 *= denom2;
  146.       adv = 0;
  147.  
  148.       remain = (LONG)multiplier * denom;
  149.       while (NotZero)
  150.       {
  151.             for (x = adv; x < size; x++)
  152.             {
  153.                   temp = (LONG)powers[x] + remain;
  154.                   powers[x] = (SHORT)(temp / denom2);
  155.                   remain = (temp - (denom2 * (LONG)powers[x])) * BASE;
  156.             }
  157.  
  158.             remain = 0;
  159.             for (x = adv; x < size; x++)
  160.             {
  161.                   temp = (LONG)powers[x] + remain;
  162.                   term[x] = (SHORT)(temp / divisor);
  163.                   remain = (temp - (divisor * (LONG)term[x])) * BASE;
  164.             }
  165.             remain = 0;
  166.  
  167.             if (sign > 0)
  168.             {
  169.                   LONG carry, sum;
  170.  
  171.                   carry = 0;
  172.                   for (x = size - 1; x >=0; x--)
  173.                   {
  174.                         sum = (LONG)pi[x] + (LONG)term[x] + carry;
  175.                         carry = 0;
  176.                         if (sum >= BASE)
  177.                         {
  178.                               carry = 1;
  179.                               sum -= BASE;
  180.                         }
  181.                         pi[x] = (SHORT)sum;
  182.                   }
  183.             }
  184.             else
  185.             {
  186.                   LONG borrow, sum;
  187.  
  188.                   borrow = 0;
  189.                   for (x = size - 1; x >= 0; x--)
  190.                   {
  191.                         sum = (LONG)pi[x] - (LONG)term[x] - borrow;
  192.                         borrow = 0;
  193.                         if (sum < 0)
  194.                         {
  195.                               borrow = 1;
  196.                               sum += BASE;
  197.                         }
  198.                         pi[x] = (SHORT)sum;
  199.                   }
  200.             }
  201.  
  202.             sign = -sign;
  203.             divisor += 2;
  204.             NotZero = 0;
  205.             for (x = adv; x < size; x++)
  206.             {
  207.                   if (powers[x])
  208.                   {
  209.                         NotZero = 1;
  210.                         break;
  211.                   }
  212.             }
  213.  
  214.             if (NotZero)
  215.             {
  216.                   while (powers[adv] == 0)
  217.                         adv++;
  218.             }
  219.             /* We can skip ones that are already 0 */
  220.       }
  221. }
  222.  
  223. int main(int argc, char *argv[])
  224. {
  225.       INDEXER x;
  226.       time_t T1, T2;
  227.  
  228.       if (argc != 2)
  229.       {
  230.             printf("I need to know how many digits to compute.\n");
  231.             exit(EXIT_FAILURE);
  232.       }
  233.  
  234.       size = atol(argv[1]);
  235.       if (size <= 0L)
  236.       {
  237.             printf("Invalid argument.\n");
  238.             exit(EXIT_FAILURE);
  239.       }
  240.  
  241.       size = ((size + BASEDIGITS - 1) / BASEDIGITS) + 1;
  242.  
  243.       pi = malloc(sizeof(SHORT) * size);
  244.       powers = malloc(sizeof(SHORT) * size);
  245.       term = malloc(sizeof(SHORT) * size);
  246.  
  247.       if ((pi == NULL) || (powers == NULL) || (term == NULL))
  248.       {
  249.             printf("Unable to allocate enough memory.\n");
  250.             exit(EXIT_FAILURE);
  251.       }
  252.  
  253.       for (x = 0; x < size; x++)
  254.             pi[x] = 0;
  255.  
  256.       T1 = time(NULL);
  257.  
  258. #if defined ARC3
  259.       arctan(8, 3, 1);
  260.       arctan(4, 7, 1);
  261. #elif defined ARC5
  262.       arctan(16, 5, 1);
  263.       arctan(4, 70, -1);
  264.       arctan(4, 99, 1);
  265. #elif defined  ARC4
  266.       arctan(12, 4, 1);
  267.       arctan(4, 20, 1);
  268.       arctan(4, 1985, 1);
  269. #elif defined  ARC10
  270.       arctan(32, 10, 1);
  271.       arctan(4, 239, -1);
  272.       arctan(16, 515, -1);
  273. #else
  274.       /* Machin formula */
  275.       arctan(16, 5, 1);
  276.       arctan(4, 239, -1);
  277. #endif
  278.  
  279.       T2 = time(NULL);
  280.  
  281.       Print(pi);
  282.  
  283. #ifdef SHOWTIME
  284.       printf("\nCalculation time %0.0lf\n", difftime(T2, T1));
  285. #endif
  286.  
  287.       return EXIT_SUCCESS;
  288. }
  289.