Zoom in and out with the mouse wheel.
[ipdf/code.git] / 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> // for PI
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
35         sum = (f(xmin) + T(4.0)*odd + T(2.0)*even + f(xmin+T(intervals)*dx))*dx / T(3.0);
36         return sum;
37
38 }
39
40 #define MAX_INTERVALS 1e9
41
42 /**
43  * main
44  */
45 int main(int argc, char ** argv)
46 {
47         printf("# intervals\terror_float\tclock_float\terror_double\tclock_double\terror_real\tclock_real\n");
48         for (uint64_t intervals = 1; intervals < (uint64_t)(MAX_INTERVALS); intervals*=10)
49         {
50                 clock_t start = clock();
51                 float error_float = Integrate<float>(f<float>, 0.0, 1.0, intervals) - (float)(M_PI);
52                 clock_t clock_float = clock() - start;
53                 start = clock();
54                 double error_double = Integrate<double>(f<double>, 0.0, 1.0, intervals) - (double)(M_PI);
55                 clock_t clock_double = clock() - start;
56                 start = clock();
57                 Real error_real = Integrate<Real>(f<Real>, 0.0, 1.0, intervals) - Real(M_PI);
58                 clock_t clock_real = clock() - start;
59         
60                 printf("%lu\t%.30f\t%li\t%.30f\t%li\t%.30f\t%li\n", intervals, error_float, clock_float, error_double, clock_double, Float(error_real), clock_real);
61         }
62 }
63

UCC git Repository :: git.ucc.asn.au