src/tests/calculatepi.cpp

   1 /**
   2  * Tester
   3  * Calculates PI using integration
   4  * Compares results obtained with float, double, and Real.
   5  */
   6 #include "main.h"
   7 #include <cmath>
   8 #include <ctime> // for performance measurements
   9
  10 /** Function to integrate
  11  *
  12  */
  13 template <class T> T f(const T & x)
  14 {
  15         // Use reference because Real might get big
  16         return T(4.0) / (T(1.0) + x*x);
  17 }
  18
  19 /**
  20  * Integrate f using the simpson rule
  21  */
  22 template <class T> T Integrate(T(*f)(const T & ), const T & xmin, const T & xmax, uint64_t intervals)
  23 {
  24         T sum = 0.0;
  25         T dx = (xmax - xmin) / T(intervals);
  26
  27         T odd = 0.0;
  28         T even = 0.0;
  29         uint64_t i = 0;
  30         for (i = 1; i < intervals-1; i+=2)
  31                 odd += f(xmin + dx*T(i));
  32         for (i = 2; i < intervals-1; i+=2)
  33                 even += f(xmin + dx*T(i));
  34         sum = (f(xmin) + T(4.0)*odd + T(2.0)*even + f(xmin+T(intervals)*dx))*dx / T(3.0);
  35         return sum;
  36
  37 }
  38
  39 #define MAX_INTERVALS 1e10
  40 /**
  41  * main
  42  */
  43 int main(int argc, char ** argv)
  44 {
  45         setbuf(stdout, NULL);
  46         setbuf(stderr, NULL);
  47         long double PI = acosl(-1.0L);
  48         printf("# intervals\terror_float\tclock_float\terror_double\tclock_double\terror_long\tclock_long");
  49         #if REAL > REAL_LONG_DOUBLE
  50                 printf("\terror_real\tclock_real\n");
  51         #else
  52                 printf("\n");
  53         #endif //REAL
  54         for (uint64_t intervals = 1; intervals < (uint64_t)(MAX_INTERVALS); intervals*=5)
  55         {
  56                 clock_t start = clock();
  57                 float error_float = Integrate<float>(f<float>, 0.0f, 1.0f, intervals) - float(PI);
  58                 clock_t clock_float = clock() - start;
  59                 start = clock();
  60                 double error_double = Integrate<double>(f<double>, 0.0, 1.0, intervals) - double(PI);
  61                 clock_t clock_double = clock() - start;
  62                 start = clock();
  63                 long double error_long = Integrate<long double>(f<long double>,0.0L,1.0L, intervals) - PI;
  64                 clock_t clock_long = clock() - start;
  65
  66                 printf("%lu\t%.30f\t%li\t%.30lf\t%li\t%.30llf\t%li", intervals, error_float, clock_float, error_double, clock_double, error_long, clock_long);
  67
  68                 #if REAL > REAL_LONG_DOUBLE
  69                         Real error_real = Integrate<Real>(f<Real>,Real(0.0L), Real(1.0L), intervals) - Real(PI);
  70                         clock_t clock_real = clock() - start;
  71                         printf("\t%.30lf\t%li\n", Float(error_real), clock_real);
  72                 #else
  73                         printf("\n");
  74                 #endif
  75         }
  76 }
  77