Inner step of multiplication in x64 assembly
[ipdf/code.git] / src / tests / qtest.cpp
1 #include <stdlib.h>
2 #include <math.h>
3 #include <stdio.h>
4 #include "main.h"
5
6 #include "real.h"
7
8 using namespace IPDF;
9
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  */
14
15
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 }
33
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));
49
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);
54         
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 }
61
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         }
70
71         size_t n = 12;
72         Real r[n];
73         r[0] = (1 << 12) + 2.0;
74         r[1] = (1 << 12) + 2.25;
75         r[2] = (16 << 8) + 1.0 + 1.0/(16<<4);
76         r[3] = (1<<24) + 2.0;
77         r[4] = (1<<24) + 3.0;
78         r[5] = 94906267.0;
79         r[6] = 94906267 + 0.25;
80         r[7] = (1<<28) - 5.5;
81         r[8] = (1<<28) - 4.5;
82         r[9] = (1<<28) + 2.0;
83         r[10] = (1<<28) + 2.25;
84         r[11] = (16<<24) + 1 + 1.0/(16<<20);
85 //      r[12] = (1<<32) + 2.0;
86 //      r[13] = (1<<32) + 2.25;
87
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