From e6e138e8baad2dda39cecbb35557710f2a550398 Mon Sep 17 00:00:00 2001
From: Sam Moore <matches@ucc.asn.au>
Date: Thu, 1 May 2014 18:55:28 +0800
Subject: [PATCH] Tester from some lecture notes by Kahan

Supposedly this is a benchmark for determining the worst accuracy of a CPU using roots to a quadratic
I think something wierd is going on; a result of infinity means there is no rounding errors but I'm not sure what results of -nan mean
and that's what I get.

The article is referenced as kahan1996ieee754 in the literature notes
---
 src/tests/qtest.cpp | 97 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 97 insertions(+)
 create mode 100644 src/tests/qtest.cpp

diff --git a/src/tests/qtest.cpp b/src/tests/qtest.cpp
new file mode 100644
index 0000000..00e309b
--- /dev/null
+++ b/src/tests/qtest.cpp
@@ -0,0 +1,97 @@
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include "main.h"
+
+#include "real.h"
+
+using namespace IPDF;
+
+/**
+ * This is a C version of Kahan's proposed benchmark for determining the worst accuracy of a CPU
+ * Taken from kahan1996iee754 - Lecture notes on the status of IEEE 754
+ */
+
+
+/**
+ * Quadradic of the form p*x^2 - 2 q*x + r == 0
+ * This function computes the roots x1 and x2 as accuractely as they are determined by {p, q, r}
+ */
+void Qdrtc(Real p, Real q, Real r, Real * x1, Real * x2)
+{
+	Real s = sqrt(q*q - p*r);
+	Real S = q + copysignf(s, q);
+	if (S == 0.f)
+	{
+		*x1 = r/p;
+		*x2 = r/p;
+		return;
+	}
+	*x1 = r/S;
+	*x2 = S/p;	
+}
+
+/**
+ * Kahan's magic Qtrial test thing
+ * Calculate roots using Qdrtc and compare to known results
+ * Return the smallest (?) of the errors
+ */
+Real Qtrial(Real r)
+{
+	Real p = r-Real(2);
+	Real q = r-Real(1);
+	Real result = Real(0);
+	Debug("Qdrtc for r = %.30lf", Float(r));
+	if (p <= Real(0))
+		Fatal("Expect r > 2");
+	else if (!((r-q)== Real(1) && (q-p) == Real(1)))
+		Fatal("r too big for Qtrial %.30lf", Float(r));
+
+	Real x1; Real x2;
+	Qdrtc(p, q, r, &x1, &x2);
+	Real e1 = -log2f(x1 - Real(1));
+	Real e2 = -log2f((x2 - Real(1)) - Real(2)/p);
+	
+	result = min(e1, e2);
+	Debug("gets %.30lf and %.30lf sig bits", Float(e1), Float(e2));
+	if (x1 < Real(1))
+		Debug(" and root %.30lf isn't at least 1", Float(x1));
+	return result;
+}
+
+int main(int argc, char ** argv)
+{
+	Real x1; Real x2;
+	Qdrtc(2,5,12,&x1,&x2);
+	if (x1 != Real(2) || x2 != Real(3))
+	{
+		Fatal("Qdrtc(2, 5, 12) failed sanity check; x1 %.30lf x2 %.30lf", Float(x1), Float(x2));
+	}
+
+	size_t n = 12;
+	Real r[n];
+	r[0] = (1 << 12) + 2.0;
+	r[1] = (1 << 12) + 2.25;
+	r[2] = (16 << 8) + 1.0 + 1.0/(16<<4);
+	r[3] = (1<<24) + 2.0;
+	r[4] = (1<<24) + 3.0;
+	r[5] = 94906267.0;
+	r[6] = 94906267 + 0.25;
+	r[7] = (1<<28) - 5.5;
+	r[8] = (1<<28) - 4.5;
+	r[9] = (1<<28) + 2.0;
+	r[10] = (1<<28) + 2.25;
+	r[11] = (16<<24) + 1 + 1.0/(16<<20);
+//	r[12] = (1<<32) + 2.0;
+//	r[13] = (1<<32) + 2.25;
+
+	Real e = +INFINITY;
+	Real t;
+	for (size_t j = 0; j < n; ++j)
+	{
+		t = Qtrial(r[j]);
+		if (t < e || t != t) e = t;
+	}
+	Debug("Worst accuracy is %.30lf sig. bits", Float(e));
+	return 0;
+}
-- 
2.20.1