1 #include <stdlib.h>
2 #include <math.h>
3 #include <stdio.h>
4 #include "main.h"
6 #include "real.h"
8 using namespace IPDF;
10 /**
11  * This is a C version of Kahan's proposed benchmark for determining the worst accuracy of a CPU
12  * Taken from kahan1996iee754 - Lecture notes on the status of IEEE 754
13  */
16 /**
17  * Quadradic of the form p*x^2 - 2 q*x + r == 0
18  * This function computes the roots x1 and x2 as accuractely as they are determined by {p, q, r}
19  */
20 void Qdrtc(Real p, Real q, Real r, Real * x1, Real * x2)
21 {
22         Real s = sqrt(q*q - p*r);
23         Real S = q + copysignf(s, q);
24         if (S == 0.f)
25         {
26                 *x1 = r/p;
27                 *x2 = r/p;
28                 return;
29         }
30         *x1 = r/S;
31         *x2 = S/p;
32 }
34 /**
35  * Kahan's magic Qtrial test thing
36  * Calculate roots using Qdrtc and compare to known results
37  * Return the smallest (?) of the errors
38  */
39 Real Qtrial(Real r)
40 {
41         Real p = r-Real(2);
42         Real q = r-Real(1);
43         Real result = Real(0);
44         Debug("Qdrtc for r = %.30lf", Float(r));
45         if (p <= Real(0))
46                 Fatal("Expect r > 2");
47         else if (!((r-q)== Real(1) && (q-p) == Real(1)))
48                 Fatal("r too big for Qtrial %.30lf", Float(r));
50         Real x1; Real x2;
51         Qdrtc(p, q, r, &x1, &x2);
52         Real e1 = -log2f(x1 - Real(1));
53         Real e2 = -log2f((x2 - Real(1)) - Real(2)/p);
55         result = min(e1, e2);
56         Debug("gets %.30lf and %.30lf sig bits", Float(e1), Float(e2));
57         if (x1 < Real(1))
58                 Debug(" and root %.30lf isn't at least 1", Float(x1));
59         return result;
60 }
62 int main(int argc, char ** argv)
63 {
64         Real x1; Real x2;
65         Qdrtc(2,5,12,&x1,&x2);
66         if (x1 != Real(2) || x2 != Real(3))
67         {
68                 Fatal("Qdrtc(2, 5, 12) failed sanity check; x1 %.30lf x2 %.30lf", Float(x1), Float(x2));
69         }
71         size_t n = 12;
72         Real r[n];
73         r = (1 << 12) + 2.0;
74         r = (1 << 12) + 2.25;
75         r = (16 << 8) + 1.0 + 1.0/(16<<4);
76         r = (1<<24) + 2.0;
77         r = (1<<24) + 3.0;
78         r = 94906267.0;
79         r = 94906267 + 0.25;
80         r = (1<<28) - 5.5;
81         r = (1<<28) - 4.5;
82         r = (1<<28) + 2.0;
83         r = (1<<28) + 2.25;
84         r = (16<<24) + 1 + 1.0/(16<<20);
85 //      r = (1<<32) + 2.0;
86 //      r = (1<<32) + 2.25;
88         Real e = +INFINITY;
89         Real t;
90         for (size_t j = 0; j < n; ++j)
91         {
92                 t = Qtrial(r[j]);
93                 if (t < e || t != t) e = t;
94         }
95         Debug("Worst accuracy is %.30lf sig. bits", Float(e));
96         return 0;
97 } UCC git Repository :: git.ucc.asn.au