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 }