"Fossies" - the Fresh Open Source Software Archive

Member "mapm_4.9.5a/PI_DEMO/pi_4.cpp" (21 Feb 2010, 4158 Bytes) of package /linux/misc/old/mapm-4.9.5a.tar.gz:

    1 
    2 /****************************************************************************/
    3 
    4 /*
    5  *      PI_4.CPP
    6  *
    7  *      DEMO 'C++' PROGRAM TO COMPUTE PI USING
    8  *      J.BORWEIN-P.BORWEIN'S ALGORITHM (1995)
    9  *
   10  */
   11 
   12 #include <stdio.h>
   13 #include <stdlib.h>
   14 #include <string.h>
   15 #include "m_apm.h"
   16 
   17 #define  FALSE  0
   18 #define  TRUE   1
   19 
   20 
   21 extern   MAPM   compute_PI_2(int);
   22 
   23 
   24 int main(int argc, char *argv[])
   25 {
   26 int      ii, dplaces;
   27 char     *buffer;
   28 MAPM     pi_mapm;            /* declare the MAPM variable ... */
   29 
   30 if (argc < 2)
   31   {
   32    fprintf(stdout,"\nUsage: pi_4  digits\t\t\t\t[Version 1.1]\n");
   33    fprintf(stdout,
   34       "       Compute PI to \'digits\' decimal places of accuracy\n");
   35 
   36    exit(4);
   37   }
   38 
   39 dplaces = atoi(argv[1]);
   40 
   41 if (dplaces < 2)
   42   {
   43    fprintf(stdout,"\nUsage: pi-4  digits \n");
   44    exit(6);
   45   }
   46 
   47 fprintf(stderr, "\nstart Borwein-Borwein PI calculation \n");
   48 
   49 pi_mapm = compute_PI_2(dplaces);
   50 
   51 fprintf(stderr, "Borwein-Borwein PI calculation complete \n\n");
   52 
   53 /*
   54  *   find out how many digits so we can
   55  *   convert the entire value into a string.
   56  *   (plus some pad for decimal point, exponent, etc)
   57  */
   58 
   59 ii = pi_mapm.significant_digits() + 12;
   60 
   61 if ((buffer = (char *)malloc(ii)) == NULL)
   62   {
   63    fprintf(stdout,"PI demo: Out of memory \n");
   64    exit(8);
   65   }
   66 
   67 
   68 /* now convert the MAPM number into a string */
   69 
   70 pi_mapm.toString(buffer, dplaces);
   71 printf("Borwein PI = [%s] \n",buffer);
   72 
   73 free(buffer);
   74 
   75 exit(0);
   76 }
   77 
   78 /****************************************************************************/
   79 
   80 /*
   81  *      Calculate PI using the J. Borwein - P. Borwein's Algorithm (1995)
   82  *
   83  *
   84  *      Init :  A0 = 6 - 4 * sqrt(2)
   85  *              B0 = sqrt(2) - 1
   86  *
   87  *      Iterate:
   88  *
   89  *                                  4
   90  *                  1 - quart( 1 - B  ) 
   91  *                                  k  
   92  *      B      =   --------------------
   93  *       k+1                        4
   94  *                  1 + quart( 1 - B  )       
   95  *                                  k
   96  *
   97  *
   98  *                                4    2k+3                   2
   99  *      A      =   A  ( 1 + B    )  - 2     B   ( 1 + B    + B    )    
  100  *       k+1        k        k+1             k+1       k+1    k+1
  101  *
  102  *
  103  *
  104  *
  105  *                  where quart(x) = sqrt(sqrt(x)) = 4th root of x
  106  *      
  107  *                     1
  108  *      A      --->  -----
  109  *       k+1          PI
  110  *
  111  */
  112 
  113 MAPM    compute_PI_2(int places)
  114 {
  115 MAPM    t1, t2, c2, a0, a1, b0, b1;
  116 int     kk, nn, dplaces, dflag;
  117 
  118 
  119 dplaces = places + 16;               // compute a few extra digits 
  120                      // in the intermediate math 
  121 
  122 m_apm_cpp_precision(dplaces);        // tell C++ our precision requirements
  123 
  124 dflag = FALSE;
  125 kk    = 0;
  126 
  127 t1 = sqrt("2.0");                    // make sure we use the MAPM sqrt
  128 a0 = 6 - 4 * t1;
  129 b0 = t1 - 1;
  130 c2 = 2;
  131 
  132 while (TRUE)
  133   {
  134    t1 = b0.ipow(4);                  // t1 = b0 ^ 4
  135    t1 = sqrt(sqrt(1 - t1));          // t1 = root_4(1 - t1)
  136    t1 = t1.round(dplaces);
  137 
  138    b1 = (1 - t1) / (1 + t1);
  139    b1 = b1.round(dplaces);
  140 
  141    t2 = 1 + b1;
  142    t2 = t2.ipow(4);                  // t2 = (1 + b1) ^ 4
  143    t2 = t2.round(dplaces);
  144 
  145    t1 = c2.ipow(2 * kk + 3);         // t1 = 2 ^ (2*kk+3)
  146 
  147    a1 = a0 * t2 - t1 * b1 * (1 + b1 + b1 * b1);
  148    a1 = a1.round(dplaces);
  149 
  150    if (dflag)
  151      break;
  152    
  153    //  compute the difference from this 
  154    //  iteration to the last one.
  155 
  156    t1 = a1 - a0;
  157 
  158    if (kk >= 3)
  159      {
  160       //  if diff == 0, we're done.
  161 
  162       if (t1 == 0)
  163         break;
  164      }
  165 
  166    /*
  167     *   if the exponent of the error term (small error like 2.47...E-65)
  168     *   is small enough, break out after the *next* a1 is calculated.
  169     *
  170     *   get the exponent of the error term, which will be negative.
  171     */
  172 
  173    nn = -(t1.exponent());
  174 
  175    //  normally, this wouldn't be here. it's nice in the demo though.
  176 
  177    fprintf(stderr, "PI now known to %d digits of accuracy ... \n",(4 * nn));
  178 
  179    if ((15 * nn) >= dplaces)
  180      dflag = TRUE;
  181 
  182    //  set up for the next iteration
  183 
  184    b0 = b1;
  185    a0 = a1;
  186    kk++;
  187   }
  188 
  189 //  round to the accuracy asked for after taking the reciprocal
  190 
  191 t1 = 1 / a1;
  192 return(t1.round(places));
  193 }
  194 
  195 /****************************************************************************/
  196