home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / snip9707.zip / PI_AGM.C < prev    next >
C/C++ Source or Header  |  1997-07-05  |  18KB  |  673 lines
  1. /* +++Date last modified: 05-Jul-1997 */
  2.  
  3. /*
  4. ** This method implements the Salamin / Brent / Gauss Arithmetic-
  5. ** Geometric Mean pi formula.
  6. **
  7. ** Let A[0] = 1, B[0] = 1/Sqrt(2)
  8. **
  9. ** Then iterate from 1 to 'n'.
  10. ** A[n] = (A[n-1] + B[n-1])/2
  11. ** B[n] = Sqrt(A[n-1]*B[n-1])
  12. ** C[n]^2 = A[n]^2 - B[n]^2  (or) C[n] = (A[n-1]-B[n-1])/2
  13. **                        n
  14. ** PI[n] = 4A[n+1]^2 / (1-(Sum (2^(j+1))*C[j]^2))
  15. **                       j = 1
  16. **
  17. ** There is an actual error calculation, but it comes out  to  slightly
  18. ** more  than double on each iteration.  I think it results in about 17
  19. ** million  correct  digits,  instead  of  16  million  if  it actually
  20. ** doubled.  PI16 generates 178,000 digits.  PI19 to  over  a  million.
  21. ** PI22 is 10 million, and PI26 to 200 million.
  22. **
  23. ** For what little it's worth, this program is placed into the public
  24. ** domain by its author, Carey Bloodworth, on September 21, 1996.
  25. **
  26. ** The math routines in this program are not general purpose routines.
  27. ** They have been optimized and written for this specific use.  The
  28. ** concepts are of course portable, but not the implementation.
  29. **
  30. ** The program run time is about O(n^1.7).  That's fairly good growth,
  31. ** compared to things such as the classic arctangents.  But not as good
  32. ** as it should be, if I used a FFT based multiplication.  Also, the 'O'
  33. ** is fairly large.  In fact, I'd guess that this program could compute
  34. ** one million digits of pi in about the same time as my previously
  35. ** posted arctangent method, in spite of this one having less than n^2
  36. ** growth.
  37. **
  38. ** The program has not been cleaned up.  It's written rather crudely
  39. ** and dirty.  The RSqrt() in particular is rather dirty, having gone
  40. ** through numerous changes and attempts at speeding it up.
  41. ** But I'm not planning on doing any more with it.  The growth isn't as
  42. ** low as I'd hoped, and until I find a faster multiplication, the
  43. ** method isn't any better than simpler arctangents.
  44. **
  45. ** I currently use a base of 10,000 simply because it made debugging
  46. ** easier.  A base of 65,536 and modifying the FastMul() to handle sizes
  47. ** that aren't powers of two would allow a bit better efficiency.
  48. */
  49.  
  50. #include <stdio.h>
  51. #include <stdlib.h>
  52. #include <math.h>
  53. #include <time.h>
  54. #include "snipmath.h"
  55.  
  56. #ifdef __TURBOC__
  57. typedef short int Indexer;
  58. #else
  59. typedef long int Indexer;
  60. #endif
  61.  
  62. typedef short int Short;
  63. typedef long  int Long;
  64.  
  65. Short *a, *b, *c, *Work1, *MulWork, *Work2, *Work3;
  66. int SIZE;
  67.  
  68. /*
  69. ** Work1 is explicitly used in Reciprocal() and RSqrt()
  70. ** Work1 is implicitly used in Divide() and Sqrt()
  71. ** Work2 is explicitly used in Divide() and Sqrt()
  72. ** Work3 is used only in the AGM and holds the previous reciprocal
  73. **  square root, allowing us to save some time in RSqrt()
  74. */
  75.  
  76.  
  77. void Copy(Short *a, Short *b, Indexer Len)
  78. {
  79.       while (Len--) a[Len] = b[Len];
  80. }
  81.  
  82. /*
  83. ** This rounds and copies a 2x mul result into a normal result
  84. ** Our number format will never have more than one unit of integer,
  85. ** and after a mul, we have two, so we need to fix that.
  86. */
  87.  
  88. void Round2x(Short *a, Short *b, Indexer Len)
  89. {
  90.       Indexer x;
  91.       Short carry;
  92.  
  93.       carry = 0;
  94.       if (b[Len+1] >= 5000)
  95.             carry = 1;
  96.       for (x = Len; x > 0; x--)
  97.       {
  98.             carry += b[x];
  99.             a[x-1] = carry % 10000;
  100.             carry /= 10000;
  101.       }
  102. }
  103.  
  104. void DivBy2(Short *n, Indexer Len)
  105. {
  106.       Indexer x;
  107.       Long temp;
  108.  
  109.       temp = 0;
  110.       for (x = 0; x < Len; x++)
  111.       {
  112.             temp = (Long)n[x]+temp*10000;
  113.             n[x] = (Short)(temp/2);
  114.             temp = temp%2;
  115.       }
  116. }
  117.  
  118. void DoCarries(Short *limit, Short *cur, Short carry)
  119. {
  120.       Long temp;
  121.  
  122.       while ((cur >= limit) && (carry != 0))
  123.       {
  124.             temp = *cur+carry;
  125.             carry = 0;
  126.             if (temp >= 10000)
  127.             {
  128.                   carry = 1;
  129.                   temp -= 10000;
  130.             }
  131.             *cur = temp;
  132.             --cur;
  133.       }
  134. }
  135.  
  136. void DoBorrows(Short *limit, Short *cur, Short borrow)
  137. {
  138.       Long temp;
  139.       while ((cur >= limit) && (borrow != 0))
  140.       {
  141.             temp = *cur-borrow;
  142.             borrow = 0;
  143.             if (temp < 0)
  144.         {
  145.             borrow = 1;
  146.             temp += 10000;
  147.         }
  148.             *cur = temp;
  149.             --cur;
  150.       };
  151. }
  152.  
  153. void PrintShort2(char *str, Short *num, Indexer Len)
  154. {
  155.       Indexer x;
  156.  
  157.       printf("%s ", str);
  158.       printf("%u.", num[0]);
  159.       for (x = 1; x < Len; x++)
  160.             printf("%04u", num[x]);
  161.       printf("\n");
  162. }
  163.  
  164. void PrintShort(char *str, Short *num, Indexer Len)
  165. {
  166.       Indexer x;
  167.       int printed = 0;
  168.  
  169.       printf("%s ", str);
  170.  
  171.       printf("%u.\n", num[0]);
  172.  
  173.       for (x = 1; x < Len; x++)
  174.       {
  175.             printf("%02d", num[x]/100);
  176.             printed += 2;
  177.             if ((printed % 1000) == 0)
  178.             {
  179.                   printf("\n\n\n\n");
  180.                   printed = 0;
  181.             }
  182.             else if ((printed % 50) == 0)
  183.             {
  184.                   printf("\n");
  185.             }
  186.             else if ((printed % 10) == 0)
  187.             {
  188.                   printf(" ");
  189.             }
  190.             printf("%02d", num[x] % 100);
  191.             printed += 2;
  192.             if ((printed % 1000) == 0)
  193.             {
  194.                   printf("\n\n\n\n");
  195.                   printed = 0;
  196.             }
  197.             else if ((printed % 50) == 0)
  198.             {
  199.                   printf("\n");
  200.             }
  201.             else if ((printed % 10) == 0)
  202.             {
  203.                   printf(" ");
  204.             }
  205.       }
  206.       printf("\n");
  207.  
  208. }
  209.  
  210. /* sum = a + b */
  211.  
  212. Short Add(Short *sum, Short *a, Short *b, Indexer Len)
  213. {
  214.       Long s, carry;
  215.       Indexer x;
  216.  
  217.       carry = 0;
  218.       for (x = Len - 1; x >= 0; x--)
  219.       {
  220.             s = (Long)a[x] + (Long)b[x] + carry;
  221.             sum[x] = (Short)(s%10000);
  222.             carry = s/10000;
  223.       }
  224.       return carry;
  225. }
  226.  
  227. /* dif = a-b */
  228.  
  229. Short Sub(Short *dif, Short *a, Short *b, Indexer Len)
  230. {
  231.       Long d, borrow;
  232.       Indexer x;
  233.  
  234.       borrow = 0;
  235.       for (x = Len - 1; x >= 0; x--)
  236.       {
  237.             d = (Long)a[x] - (Long)b[x] - borrow;
  238.             borrow = 0;
  239.             if (d < 0)
  240.             {
  241.                   borrow = 1;
  242.                   d += 10000;
  243.             }
  244.             dif[x] = (Short)d;
  245.       }
  246.       return borrow;
  247. }
  248.  
  249. void Negate(Short *num, Indexer Len)
  250. {
  251.       Indexer x;Long d, borrow;
  252.  
  253.       borrow = 0;
  254.       for (x = Len - 1; x >= 0; x--)
  255.       {
  256.             d = 0 - num[x] - borrow;
  257.             borrow = 0;
  258.             if (d < 0)
  259.             {
  260.                   borrow = 1;
  261.                   d += 10000;
  262.             }
  263.             num[x] = (Short)d;
  264.       }
  265. }
  266.  
  267. /* prod = a*b.  prod should be twice normal length */
  268.  
  269. void Mul(Short *prod, Short *a, Short *b, Indexer Len)
  270. {
  271.       Long p;
  272.       Indexer ia, ib, ip;
  273.  
  274.       if ((prod == a) || (b == prod))
  275.       {
  276.             printf("MUL product can't be one of the other arguments\n");
  277.             exit(1);
  278.       }
  279.  
  280.       for (ip = 0; ip < Len * 2; ip++)
  281.             prod[ip] = 0;
  282.       for (ib = Len - 1; ib >= 0; ib--)
  283.       {
  284.             if (b[ib] == 0)
  285.                   continue;
  286.             ip = ib + Len;
  287.             p = 0;
  288.             for (ia = Len - 1; ia >= 0; ia--)
  289.             {
  290.                   p = (Long)a[ia]*(Long)b[ib] + p + prod[ip];
  291.                   prod[ip] = p%10000;
  292.                   p = p/10000;
  293.                   ip--;
  294.             }
  295.             while ((p) && (ip >= 0))
  296.             {
  297.                   p += prod[ip];
  298.                   prod[ip] = p%10000;
  299.                   p = p/10000;
  300.                   ip--;
  301.             }
  302.       }
  303. }
  304.  
  305. /*
  306. ** This is based on the simple O(n^1.585) method, although my
  307. ** growth seems to be closer to O(n^1.7)
  308. **
  309. ** It's fairly simple.  a*b is: a2b2(B^2+B)+(a2-a1)(b1-b2)B+a1b1(B+1)
  310. **
  311. ** For a = 4711 and b = 6397, a2 = 47 a1 = 11 b2 = 63 b1 = 97 Base = 100
  312. **
  313. ** If we did that the normal way, we'd do
  314. ** a2b2 = 47*63 = 2961
  315. ** a2b1 = 47*97 = 4559
  316. ** a1b2 = 11*63 =  693
  317. ** a1b1 = 11*97 = 1067
  318. **
  319. ** 29 61
  320. **    45 59
  321. **     6 93
  322. **       10 67
  323. ** -----------
  324. ** 30 13 62 67
  325. **
  326. ** Or, we'd need N*N multiplications.
  327. **
  328. ** With the 'fractal' method, we compute:
  329. ** a2b2 = 47*63 = 2961
  330. ** (a2-b1)(b1-b2) = (47-11)(97-63) = 36*34 = 1224
  331. ** a1b1 = 11*97 = 1067
  332. **
  333. ** 29 61
  334. **    29 61
  335. **    12 24
  336. **    10 67
  337. **       10 67
  338. ** -----------
  339. ** 30 13 62 67
  340. **
  341. ** We need only 3 multiplications, plus a few additions.  And of course,
  342. ** additions are a lot simpler and faster than multiplications, so we
  343. ** end up ahead.
  344. **
  345. ** The way it is written requires that the size to be a power of two.
  346. ** That's why the program itself requires the size to be a power of two.
  347. ** There is no actual requirement in the method, it's just how I did it
  348. ** because I would be working close to powers of two anyway, and it was
  349. ** easier.
  350. */
  351.  
  352. void FastMul(Short *prod, Short *a, Short *b, Indexer Len)
  353. {
  354.       Indexer x, HLen;
  355.       int SignA, SignB;
  356.       Short *offset;
  357.       Short *NextProd;
  358.  
  359.       /*
  360.       ** This is the threshold where it's faster to do it the ordinary way
  361.       ** On my computer, it's 16.  Yours may be different.
  362.       */
  363.  
  364.       if (Len <= 16 )
  365.       {
  366.             Mul(prod, a, b, Len);
  367.             return;
  368.       }
  369.  
  370.       HLen = Len/2;
  371.       NextProd = prod + 2*Len;
  372.  
  373.       SignA = Sub(prod, a, a + HLen, HLen);
  374.       if (SignA)
  375.       {
  376.             SignA = 1;
  377.             Negate(prod, HLen);
  378.       }
  379.       SignB = Sub(prod + HLen, b + HLen, b, HLen);
  380.       if (SignB)
  381.       {
  382.             SignB = 1;
  383.             Negate(prod + HLen, HLen);
  384.       }
  385.  
  386.       FastMul(NextProd, prod, prod + HLen, HLen);
  387.       for (x = 0; x < 2 * Len; x++)
  388.             prod[x] = 0;
  389.       offset = prod + HLen;
  390.       if (SignA == SignB)
  391.             DoCarries(prod, offset - 1, Add(offset, offset, NextProd, Len));
  392.       else  DoBorrows(prod, offset - 1, Sub(offset, offset, NextProd, Len));
  393.  
  394.       FastMul(NextProd, a, b, HLen);
  395.       offset = prod + HLen;
  396.       DoCarries(prod, offset - 1, Add(offset, offset, NextProd, Len));
  397.       Add(prod, prod, NextProd, Len);
  398.  
  399.       FastMul(NextProd, a + HLen, b + HLen, HLen);
  400.       offset = prod + HLen;
  401.       DoCarries(prod, offset - 1, Add(offset, offset, NextProd, Len));
  402.       offset = prod + Len;
  403.       DoCarries(prod, offset - 1, Add(offset, offset, NextProd, Len));
  404. }
  405.  
  406. /*
  407. ** Compute the reciprocal of 'a'.
  408. ** x[k+1] = x[k]*(2-a*x[k])
  409. */
  410.  
  411. void Reciprocal(Short *r, Short *n, Indexer Len)
  412. {
  413.       Indexer x, SubLen;
  414.       int iter;
  415.       double t;
  416.  
  417.       /* Estimate our first reciprocal */
  418.  
  419.       for (x = 0; x < Len; x++)
  420.             r[x] = 0;
  421.       t = n[0] + n[1]*1.0e-4 + n[2]*1.0e-8;
  422.       if (t == 0.0)
  423.             t = 1;
  424.       t = 1.0/t;
  425.       r[0] = t;
  426.       t = (t - floor(t))*10000.0;
  427.       r[1] = t;
  428.       t = (t - floor(t))*10000.0;
  429.       r[2] = t;
  430.  
  431.       iter = log(SIZE)/log(2.0) + 1;
  432.       SubLen = 1;
  433.       while (iter--)
  434.       {
  435.             SubLen *= 2;
  436.             if (SubLen > SIZE)
  437.                   SubLen = SIZE;
  438.             FastMul(MulWork, n, r, SubLen);
  439.             Round2x(Work1, MulWork, SubLen);
  440.  
  441.             MulWork[0] = 2;
  442.             for (x = 1; x < Len; x++)
  443.                   MulWork[x] = 0;
  444.             Sub(Work1, MulWork, Work1, SubLen);
  445.  
  446.             FastMul(MulWork, r, Work1, SubLen);
  447.             Round2x(r, MulWork, SubLen);
  448.       }
  449. }
  450.  
  451. /*
  452. ** Computes the reciprocal of the square root of 'a'
  453. ** y[k+1] = y[k]*(3-a*y[k]^2)/2
  454. **
  455. **
  456. ** We can save quite a bit of time by using part of our previous
  457. ** results!  Since the number is converging onto a specific value,
  458. ** at least part of our previous result will be correct, so we
  459. ** can skip a bit of work.
  460. */
  461.  
  462. void RSqrt(Short *r, Short *n, Indexer Len, Indexer SubLen)
  463. {
  464.       Indexer x;
  465.       int iter;
  466.       double t;
  467.  
  468.       /* Estimate our first reciprocal square root */
  469.  
  470.       if (SubLen < 2 )
  471.       {
  472.             for (x = 0; x < Len; x++)
  473.                   r[x] = 0;
  474.             t = n[0] + n[1]*1.0e-4 + n[2]*1.0e-8;
  475.             if (t == 0.0)
  476.                   t = 1;
  477.             t = 1.0/sqrt(t);
  478.             r[0] = t;
  479.             t = (t - floor(t))*10000.0;
  480.             r[1] = t;
  481.             t = (t - floor(t))*10000.0;
  482.             r[2] = t;
  483.       }
  484.  
  485.       /*
  486.       ** PrintShort2("\n  ~R: ", r, SIZE);
  487.       */
  488.  
  489.       if (SubLen > SIZE)
  490.             SubLen = SIZE;
  491.       iter = SubLen;
  492.       while (iter <= SIZE*2)
  493.       {
  494.             FastMul(MulWork, r, r, SubLen);
  495.             Round2x(Work1, MulWork, SubLen);
  496.             FastMul(MulWork, n, Work1, SubLen);
  497.             Round2x(Work1, MulWork, SubLen);
  498.  
  499.             MulWork[0] = 3;
  500.             for (x = 1; x < SubLen; x++)
  501.                   MulWork[x] = 0;
  502.             Sub(Work1, MulWork, Work1, SubLen);
  503.  
  504.             FastMul(MulWork, r, Work1, SubLen);
  505.             Round2x(r, MulWork, SubLen);
  506.             DivBy2(r, SubLen);
  507.  
  508.             /*
  509.             printf("%3d", SubLen);
  510.             PrintShort2("R: ", r, SubLen);
  511.             printf("%3d", SubLen);
  512.             PrintShort2("R: ", r, Len);
  513.             */
  514.  
  515.             SubLen *= 2;
  516.             if (SubLen > SIZE)
  517.                   SubLen = SIZE;
  518.             iter *= 2;
  519.       }
  520. }
  521.  
  522. /*
  523. ** d = a/b by computing the reciprocal of b and then multiplying
  524. ** that by a.  It's faster this way because the normal digit by
  525. ** digit method has N^2 growth, and this method will have the same
  526. ** growth as our multiplication method.
  527. */
  528.  
  529. void Divide(Short *d, Short *a, Short *b, Indexer Len)
  530. {
  531.       Reciprocal(Work2, b, Len);
  532.       FastMul(MulWork, a, Work2, Len);
  533.       Round2x(d, MulWork, Len);
  534. }
  535.  
  536. /*
  537. ** r = sqrt(n) by computing the reciprocal of the square root of
  538. ** n and then multiplying that by n
  539. */
  540.  
  541. void Sqrt(Short *r, Short *n, Indexer Len, Indexer SubLen)
  542. {
  543.       RSqrt(Work2, n, Len, SubLen);
  544.       FastMul(MulWork, n, Work2, Len);
  545.       Round2x(r, MulWork, Len);
  546. }
  547.  
  548. void usage(void)
  549. {
  550.       puts("This program requires the number of digits/4 to calculate");
  551.       puts("This number must be a power of two.");
  552.       exit(-1);
  553. }
  554.  
  555. int main(int argc,char *argv[])
  556. {
  557.       Indexer x;
  558.       Indexer SubLen;
  559.       double Pow2, tempd, carryd;
  560.       int AGMIter;
  561.       int NeededIter;
  562.       clock_t T0, T1, T2, T3;
  563.  
  564.       if (argc < 2)
  565.             usage();
  566.  
  567.       SIZE = atoi(argv[1]);
  568.       if (!ispow2(SIZE))
  569.             usage();
  570.  
  571.       a = (Short *)malloc(sizeof(Short)*SIZE);
  572.       b = (Short *)malloc(sizeof(Short)*SIZE);
  573.       c = (Short *)malloc(sizeof(Short)*SIZE);
  574.       Work1 = (Short *)malloc(sizeof(Short)*SIZE);
  575.       Work2 = (Short *)malloc(sizeof(Short)*SIZE);
  576.       Work3 = (Short *)malloc(sizeof(Short)*SIZE);
  577.       MulWork = (Short *)malloc(sizeof(Short)*SIZE*4);
  578.  
  579.       if ((a ==NULL) || (b == NULL) || (c == NULL) || (Work1 == NULL) ||
  580.           (Work2 == NULL) || (MulWork == NULL))
  581.       {
  582.             printf("Unable to allocate memory\n");exit(1);
  583.       }
  584.  
  585.       NeededIter = log(SIZE)/log(2.0) + 1;
  586.       Pow2 = 4.0;
  587.       AGMIter = NeededIter + 1;
  588.       SubLen = 1;
  589.  
  590.       T0 = clock();
  591.  
  592.       for (x = 0; x < SIZE; x++)
  593.             a[x] = b[x] = c[x] = Work3[x] = 0;
  594.       a[0] = 1;
  595.       Work3[0] = 2;
  596.       RSqrt(b, Work3, SIZE, 1);
  597.       c[0] = 1;
  598.  
  599.       T1 = clock();
  600.  
  601.       /*
  602.       printf("AGMIter %d\n", AGMIter);
  603.       PrintShort2("a AGM: ", a, SIZE);
  604.       PrintShort2("b AGM: ", b, SIZE);
  605.       PrintShort2("C sum: ", c, SIZE);
  606.       */
  607.  
  608.       while (AGMIter--)
  609.       {
  610.             Sub(Work1, a, b, SIZE);                /* w = (a-b)/2      */
  611.             DivBy2(Work1, SIZE);
  612.             FastMul(MulWork, Work1, Work1, SIZE);  /* m = w*w          */
  613.             Round2x(Work1, MulWork, SIZE);
  614.  
  615.             carryd = 0.0;                         /* m = m* w^(J+1)   */
  616.             for (x = SIZE - 1; x >= 0; x--)
  617.             {
  618.                   tempd = Pow2*Work1[x] + carryd;
  619.                   carryd = floor(tempd / 10000.0);
  620.                   Work1[x] = tempd - carryd*10000.0;
  621.             }
  622.             Pow2 *= 2.0;
  623.             Sub(c, c, Work1, SIZE);                /* c = c - m        */
  624.  
  625.             /* Save some time on last iter */
  626.  
  627.             if (AGMIter != 0)
  628.                   FastMul(MulWork, a, b, SIZE);    /* m = a*b          */
  629.             Add(a, a, b, SIZE);                    /* a = (a+b)/2      */
  630.             DivBy2(a, SIZE);
  631.             if (AGMIter != 0)
  632.             {
  633.                   Round2x(Work3, MulWork, SIZE);
  634.                   Sqrt(b, Work3, SIZE, SubLen); /* b = sqrt(a*b) */
  635.             }
  636.             /*
  637.             printf("AGMIter %d\n", AGMIter);
  638.             PrintShort2("a AGM: ", a, SIZE);
  639.             PrintShort2("b AGM: ", b, SIZE);
  640.             PrintShort2("C sum: ", c, SIZE);
  641.             */
  642.             SubLen *= 2;if (SubLen > SIZE) SubLen = SIZE;
  643.       }
  644.  
  645.       T2 = clock();
  646.  
  647.       FastMul(MulWork, a, a, SIZE);
  648.       Round2x(a, MulWork, SIZE);
  649.       carryd = 0.0;
  650.       for (x = SIZE - 1; x >= 0; x--)
  651.       {
  652.             tempd = 4.0*a[x] + carryd;
  653.             carryd = floor(tempd / 10000.0);
  654.             a[x] = tempd - carryd*10000.0;
  655.       }
  656.  
  657.       Divide(b, a, c, SIZE);
  658.       T3 = clock();
  659.  
  660.       PrintShort("\nPI = ", b, SIZE);
  661.  
  662.       printf("\nTotal Execution time: %12.4lf\n",
  663.          (double)(T3 - T0) / CLOCKS_PER_SEC);
  664.       printf("Setup time          : %12.4lf\n",
  665.          (double)(T1 - T0) / CLOCKS_PER_SEC);
  666.       printf("AGM Calculation time: %12.4lf\n",
  667.          (double)(T2 - T1) / CLOCKS_PER_SEC);
  668.       printf("Post calc time      : %12.4lf\n",
  669.          (double)(T3 - T2) / CLOCKS_PER_SEC);
  670.  
  671.       return 0;
  672. }
  673.