home *** CD-ROM | disk | FTP | other *** search
/ Night Owl 25 / nopv25.iso / 040A / SNIP9707.ZIP / PI.C < prev    next >
Encoding:
C/C++ Source or Header  |  1997-07-05  |  4.2 KB  |  165 lines
  1. /* +++Date last modified: 05-Jul-1997 */
  2.  
  3. /*
  4. **  PI.C - Computes Pi to an arbitrary number of digits
  5. **
  6. **  Uses far arrays when/where required so may be compiled in any memory model
  7. **
  8. **  The formula that most use (including the one in the Snippets) is the
  9. **  classic Machin formula, with the Gregory series.
  10. **  
  11. **  pi=16arctan(1/5)-4arctan(1/239)
  12. **  
  13. **  And use the traditional Gregory arctangent series to calculate the
  14. **  arctangents. That's the:
  15. **  
  16. **  arctan(x)=x-(x^3)/3+(x^5)/5-(x^7)/7.....
  17. **  
  18. **  With 1/5 and 1/239, it would be:
  19. **  
  20. **  arctan(x)=1/5-1/(3*5^3)+1/(5*5^5)-1/(7*5^7)...
  21. **  arctan(x)=1/239-1/(3*239^3)+1/(5*239^5)-1/(7*239^7)....
  22. **  
  23. **  Doing the multi-precision isn't too hard, since we don't really need to
  24. **  have a general purpose math package.  We can hardwire it all.
  25. **
  26. **  Due to the % operator, ms[i] < (temp * (239**2)) so
  27. **  temp < 3759 and i < 1879, it fails at the 1879th term which translates to
  28. **  1879 * log10(239**2) == 8941th decimal.
  29. **
  30. **  In practice we get a few more digits, (2 -> 8943th)
  31. */
  32.  
  33. #include <stdio.h>
  34. #include <stdlib.h>
  35. #include "big_mall.h"
  36.  
  37. long kf, ks;
  38. long FAR *mf, FAR *ms;
  39. long cnt, n, temp, nd;
  40. long i, line;
  41. long col, col1;
  42. long loc, FAR stor[4096];
  43.  
  44. static void PASCAL shift(long FAR *l1, long FAR *l2, long lp, long lmod)
  45. {
  46.       long k;
  47.  
  48.       k = ((*l2) > 0 ? (*l2) / lmod: -(-(*l2) / lmod) - 1);
  49.       *l2 -= k * lmod;
  50.       *l1 += k * lp;
  51. }
  52.  
  53. static void PASCAL yprint(long m)
  54. {
  55.       if (cnt<n)
  56.       {
  57.             if (++col == 11)
  58.             {
  59.                   col = 1;
  60.                   if (++col1 == 5)
  61.                   {
  62.                         col1 = 0;
  63.                         printf("\nL%04ld:", ++line);
  64.                         printf("%4ld",m%10);
  65.                   }
  66.                   else  printf("%3ld",m%10);
  67.             }
  68.             else  printf("%ld",m);
  69.             cnt++;
  70.       }
  71. }
  72.  
  73. static void PASCAL xprint(long m)
  74. {
  75.       long ii, wk, wk1;
  76.  
  77.       if (m < 8)
  78.       {
  79.             for (ii = 1; ii <= loc; )
  80.                   yprint(stor[(int)(ii++)]);
  81.             loc = 0;
  82.       }
  83.       else
  84.       {
  85.             if (m > 9)
  86.             {
  87.                   wk = m / 10;
  88.                   m %= 10;
  89.                   for (wk1 = loc; wk1 >= 1; wk1--)
  90.                   {
  91.                         wk += stor[(int)wk1];
  92.                         stor[(int)wk1] = wk % 10;
  93.                         wk /= 10;
  94.                   }
  95.             }
  96.       }
  97.       stor[(int)(++loc)] = m;
  98. }
  99.  
  100. static void PASCAL memerr(int errno)
  101. {
  102.         printf("\a\nOut of memory error #%d\n", errno);
  103.         if (2 == errno)
  104.                 BigFree(mf);
  105.         _exit(2);
  106. }
  107.  
  108. int main(int argc, char *argv[])
  109. {
  110.       int i=0;
  111.       char *endp;
  112.  
  113.       stor[i++] = 0;
  114.       if (argc < 2)
  115.       {
  116.             puts("\aUsage: PI <number_of_digits>");
  117.             return(1);
  118.       }
  119.       n = strtol(argv[1], &endp, 10);
  120.       mf = BigMalloc((Size_T)(n + 3L), (Size_T)sizeof(long));
  121.       if (NULL == mf)
  122.             memerr(1);
  123.       ms = BigMalloc((Size_T)(n + 3L), (Size_T)sizeof(long));
  124.       if (NULL == ms)
  125.             memerr(2);
  126.       printf("\nApproximation of PI to %ld digits\n", (long)n);
  127.       cnt = 0;
  128.       kf = 25;
  129.       ks = 57121L;
  130.       mf[1] = 1L;
  131.       for (i = 2; i <= (int)n; i += 2)
  132.       {
  133.             mf[i] = -16L;
  134.             mf[i+1] = 16L;
  135.       }
  136.       for (i = 1; i <= (int)n; i += 2)
  137.       {
  138.             ms[i] = -4L;
  139.             ms[i+1] = 4L;
  140.       }
  141.       printf("\nL%04ld: 3.", ++line);
  142.       while (cnt < n)
  143.       {
  144.             for (i = 0; ++i <= (int)n - (int)cnt; )
  145.             {
  146.                   mf[i] *= 10L;
  147.                   ms[i] *= 10L;
  148.             }
  149.             for (i =(int)(n - cnt + 1); --i >= 2; )
  150.             {
  151.                   temp = 2 * i - 1;
  152.                   shift(&mf[i - 1], &mf[i], temp - 2, temp * kf);
  153.                   shift(&ms[i - 1], &ms[i], temp - 2, temp * ks);
  154.             }
  155.             nd = 0;
  156.             shift((long FAR *)&nd, &mf[1], 1L, 5L);
  157.             shift((long FAR *)&nd, &ms[1], 1L, 239L);
  158.             xprint(nd);
  159.       }
  160.       printf("\n\nCalculations Completed!\n");
  161.       BigFree(ms);
  162.       BigFree(mf);
  163.       return(0);
  164. }
  165.