home *** CD-ROM | disk | FTP | other *** search
/ Archive Magazine CD 1995 / Archive Magazine CD 1995.iso / discs / prog_disc / volume_5 / issue_06 / benchmarks / c_source / newpi < prev    next >
Encoding:
Text File  |  1991-04-29  |  7.0 KB  |  224 lines
  1. /* C version of a program to calculate PI. Version 4,  87-11-29
  2.                  B A Wichmann.
  3.  
  4.    Special notes for the C version:
  5.  
  6.    a) Converted by hand from Standard Pascal.
  7.  
  8.    b) This is the 32-bit version.
  9.  
  10.    c) It calculates to a fixed number (500) digits, to remove
  11.       terminal i/o variability.
  12.  
  13.   This program is an integer benchmark, that is, a program which reflects
  14.   the speed of a computer performing integer calculations. The computation
  15.   used is an 'actual' one, rather than synthetic (as used by most PC
  16.   benchmarks). The program can be used to compare the processing speeds
  17.   of 16-bit and 32-bit implementations.
  18.  
  19.     (The following changes are needed to the program text for
  20.      successful execution on a 16-bit system:
  21.  
  22.                   Bits16 = true;
  23.                   Radix = 10;
  24.                   DRadix = 1;
  25.                   MaxT = 1000;
  26.                   last declaration, 251 changed to 1001
  27.      Change the output statement to add ',1', see comment.)
  28.  
  29.   Since the main computation uses integer arithmetic, the results are
  30.   repeatable (except for the last line which gives the rounding error
  31.   estimate which could vary by one digit).
  32.   The program can be timed by adding timing statements to the program
  33.   or by timing the whole program (ensuring that the input of the data
  34.   does not cause any delay).
  35.   The time taken for the 32-bit version of the program is about:
  36.       = K * Digits * ( Digits + 5.5 )
  37.   where K depends upon the machine and is approximately 20 machine
  38.   instruction times. The machine instruction time for an interpretive
  39.   system is that of a single interpreted instruction. For example,
  40.   on a BBC 'B' the interpreter time is about 200 microseconds per
  41.   instruction.  (Digits is the number of digits requested.)
  42.   Note that the program can be converted to other languages quite
  43.   easily. Such conversions should adhere to the form of the Pascal
  44.   version if comparisons are to be meaningful. (If fact, the program
  45.   was first written in BASIC.) The 16-bit version is much slower
  46.   since the radix for the computation has to be reduced from 10000
  47.   to 10. Version 4 of this program contains additional checks to
  48.   ensure correct execution since experience has shown that this
  49.   is necessary, especially if the program is translated into other
  50.   languages.
  51.  
  52. */
  53.  
  54. #include <stdio.h>
  55. #include <time2.h>
  56.  
  57.  
  58. #define   Bits16  0 /* false for 32-bit, true for 16-bit version */
  59. #define   MaxD    1000 /* Maximum number of digits permitted */
  60. #define   Radix  10000 /* Radix for calculation of 4 digits  */
  61.                   /* NB. Radix of 5 digits overflows on 32-bit integers  */
  62. #define   DRadix  4    /* Number of zeros in Radix, Radix=10^DRadix */
  63. #define   MaxT   250   /* MaxD / DRadix */
  64.  
  65. typedef long int CompInt;
  66.  
  67. static int   Digits;  /* Number of digits required */
  68. static CompInt   Ar[MaxT+1];
  69. static CompInt   Br[MaxT+1];
  70. static CompInt   Sr[MaxT+1];
  71. static char   St[MaxT+1][DRadix+1];
  72.  
  73. void Check(position, digit)
  74. int position,digit;
  75.  
  76.    {
  77.    if (position < Digits - 6)
  78.       {
  79.       if (St[(position-1) / DRadix + 2] [((position-1) % DRadix + 1)] !=
  80.          (digit+ 48))
  81.          printf("Incorrect value, position %d, digit%d\n",
  82.                                    position, digit);
  83.       }
  84.    }
  85.  
  86. main()
  87. /* This program calculates PI to arbitrary accuracy */
  88.  
  89.       /* PI= 8*arctan(1/3) + 4*arctan(1/7)
  90.           =2.4(1+2*1/(3*10)+2*4*1/(3*5*10^2)+ ...)
  91.            +.56(1+2*2/(3*100)+2*4/(3*5)*(2/100)^2+ ...)  */
  92.  
  93. {
  94.  
  95. CompInt   A, B, Sc;
  96. unsigned short int   Terms;
  97. unsigned short int   d,j;
  98. CompInt   x, Fivex, y, i;
  99. unsigned short int   k, kminus1; /* Limit is MaxT + 1 */
  100. clock_t start, end;
  101.  
  102. printf("Working... \n");
  103.  
  104. start=clock();
  105.  
  106. /*   printf("Number of digits required\n");        */
  107. /*   scanf(" %2d",&Digits);                        */
  108.    Digits = 500;
  109.    if ((Digits > 1000) || (Bits16 && (Digits > 550)))
  110.       printf("Too many digits\n");
  111.    Terms = Digits / DRadix; /* = terms to be used of Acc */
  112.    for (j = 1; j <= Terms; j++)
  113.       {
  114.       Ar[j] = 0; Br[j] = 0; Sr[j] = 0;
  115.       }
  116.    Sr[1] = 2;
  117.    Ar[1] = 2;
  118.    if (Bits16) {
  119.       if ((Radix !=10) || (DRadix !=0) || (MaxT != 1000)) {
  120.          printf("Wrong program specification\n");
  121.       }
  122.       Sr[2] = 9;
  123.       Sr[3] = 6;
  124.       Ar[2] = 4;
  125.       Br[2] = 5;
  126.       Br[3] = 6;
  127.    }
  128.    else {
  129.       Sr[2] = 9600;    /* = .96 * Radix */
  130.       Ar[2] = 4000;    /* = .4  * Radix */
  131.       Br[2] = 5600;    /* = .56 * Radix */
  132.       i = 2;  k = 1;
  133.    };
  134.    do {
  135.       A = 0;  B = 0;
  136.       if (k==1)
  137.          kminus1 = 1;
  138.       else
  139.          kminus1 = k-1;
  140.  
  141.    /* k records the number of initial zeros in Ar[] (&& hence Br[] ). */
  142.  
  143.       for (j = Terms; j >= kminus1; j--)
  144.          {
  145.          /* Multiply Ar[] && Br[] by i */
  146.          Ar[j] = (Ar[j]*i)+A;
  147.          Br[j] = (Br[j]*i)+B;
  148.          A     = Ar[j] / Radix;
  149.          Ar[j] = Ar[j]-(A*Radix);
  150.          B     = Br[j] / Radix;
  151.          Br[j] = Br[j]-(B*Radix);
  152.          } /* j*/;
  153.  
  154.       if (DRadix == 1){
  155.          x = x+1;
  156.          for (j=Terms; j <= kminus1+1; j--){
  157.              Br[j] = Br[j-1];
  158.              Ar[j] = Ar[j-1];
  159.              }
  160.          Br[kminus1] = 0;
  161.          Ar[kminus1] = 0;
  162.       }
  163.       x = 10*(i+1);
  164.       Fivex = 5*x;
  165.       A = 0;  B = 0;
  166.       for (j = kminus1; j <=Terms; j++)
  167.          {
  168.          /* Divide Ar[] by x && Br[] by Fivex */
  169.          y = (Ar[j]+A) / x;
  170.          A = Radix*((Ar[j]+A)-(y*x));
  171.          Ar[j] = y;
  172.          y = (Br[j]+B) / Fivex;
  173.          B = Radix*((Br[j]+B)-(Fivex*y));
  174.          Br[j] = y;
  175.          Sr[j] = Sr[j]+Ar[j]+Br[j];
  176.          }/*j*/;
  177.       if (Ar[k]==0){
  178.          k++;}
  179.       i += 2;
  180.    }
  181.    while  (k != Terms+1);
  182.    /* Calculation complete, perform conversion for listing */
  183.    Sc = 0;
  184.    for (j = Terms; j>=1; j--)
  185.       {
  186.       /* Normalise Sr */
  187.       Sr[j] = Sr[j]+Sc;
  188.       Sc = Sr[j] / Radix;
  189.       Sr[j] = Sr[j]-(Sc*Radix);
  190.       for (d = DRadix; d >= 1; d--)
  191.          {
  192.          St[j][d] = 48+ (Sr[j] % 10);
  193.          if ((St[j][d]<48) || (St[j][d] >57))
  194.             printf("Wrong program specification");
  195.          Sr[j] = Sr[j] / 10;
  196.         /* printf("%c",St[j][d]);*/
  197.          } /*d*/
  198.       } /*j*/
  199.    Check(  1, 1); Check( 41, 6); Check( 81, 8);
  200.    Check(121, 0); Check(161, 8); Check(201, 4);
  201.    Check(241, 2); Check(281, 1); Check(321, 4);
  202.    Check(361, 0); Check(401, 3); Check(441, 3);
  203.    Check(481, 8); Check(521, 8); Check(561, 0);
  204.    Check(601, 0); Check(641, 7); Check(681, 1);
  205.    Check(721, 8); Check(761, 4); Check(801, 5);
  206.    Check(841, 2); Check(881, 8); Check(921, 1);
  207.    Check(961, 1);
  208.    printf("   3.");
  209.    for (j = 2; j <= Terms;j++)
  210.       {
  211.       if ((Bits16 && (j % 4 == 2)) || !(Bits16))
  212.          printf(" ");
  213.       printf("%c%c%c%c",St[j][1],St[j][2],St[j][3],St[j][4]);  /* printf(St[j,1]); for 16-bit version */
  214.       if (((DRadix == 1) && (j % 40 == 1)) ||
  215.          ((DRadix == 4) && (j % 10 == 0)))
  216.          printf("\n");
  217.       }
  218.    printf("\n");
  219.    printf("Last digits inaccurate: add about %5.0f \n",
  220.        (0.78*Digits+0.0001*(Digits*Digits)) );
  221.    end=clock();
  222.    printf("The timing for this PI test, which tests integer calculation, was: %f\n", (end-start) / CLK_TCK);
  223. }
  224.