[goffice] Tests: add test program for approximating constants as doubles.

[Morten Welinder <mortenw src gnome org>]

commit 451e3ea2e9dd512fd8ed944e86fed2c4904de829
Author: Morten Welinder <terra gnome org>
Date:   Mon Mar 7 21:55:57 2016 -0500

    Tests: add test program for approximating constants as doubles.
    
    See comments in the program.

 tests/.gitignore  |    2 +
 tests/Makefile.am |    6 +-
 tests/constants.c |  202 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 208 insertions(+), 2 deletions(-)
---
diff --git a/tests/.gitignore b/tests/.gitignore
index fbb41c5..14e1d9e 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -5,6 +5,7 @@
 *.o
 .deps
 .libs
+*~
 Makefile
 Makefile.in
 go-demo
@@ -14,3 +15,4 @@ shapes-demo
 test-format
 test-math
 test-quad
+constants
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 56de845..5429508 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -1,4 +1,4 @@
-check_PROGRAMS=test-quad test-math test-format
+check_PROGRAMS=test-quad test-math test-format constants
 if WITH_GTK
 check_PROGRAMS += pie-demo go-demo shapes-demo mf-demo
 endif
@@ -9,6 +9,9 @@ AM_CFLAGS = $(GOFFICE_CFLAGS)
 
 TESTS = test-quad test-math
 
+constants_LDADD = $(GOFFICE_PLUGIN_LIBADD)
+constants_SOURCES = constants.c
+
 pie_demo_LDADD = $(GOFFICE_PLUGIN_LIBADD)
 pie_demo_SOURCES = pie-demo.c
 
@@ -31,4 +34,3 @@ test_format_LDADD = $(GOFFICE_PLUGIN_LIBADD)
 test_format_SOURCES = test-format.c
 
 EXTRA_DIST = go-demo.ui
-
diff --git a/tests/constants.c b/tests/constants.c
new file mode 100644
index 0000000..3a4dc0f
--- /dev/null
+++ b/tests/constants.c
@@ -0,0 +1,202 @@
+#include <goffice/goffice.h>
+
+// This program answers the question of whether a constant, c, is more
+// accurately represented as a double than its inverse, 1/c.
+//
+// The purpose is to help squeeze a few more bits out of calculations
+// that involve multiplying or dividing by the constant.
+//
+// Consider the density function for the normal distribution:
+//     d(x) = (1/sqrt(2*pi)) * exp(-x^2/2)
+// As written, this multiplies by a constant, 1/sqrt(2*pi), whose value
+// we could compute and state with enough decimals to guarantee that we
+// end up using the best possible value of type double.  But we could
+// also use division by the constant sqrt(2*pi).  Which is best?
+//
+// The traditional answer has been that division is expensive, so avoid
+// it, but that isn't really true anymore.  Here we are only interested
+// in accuracy.
+//
+// Assumptions:
+// 1. double is base-2 floating-point with 53 mantissa bits (including
+//    the implied leading bit).
+// 2. we want to minimize the error on the product or quotient.
+// 3. we know nothing special about the second factor that would
+//    force one choice or the other.
+//
+// For sqrt(2pi) the answer is that it is ever so slightly better to use
+// multiplication.  But it's a very close call.  But to get the right value,
+// one cannot start with the best double approximating sqrt(2pi) and compute
+// 1/sqrt(2pi) from that.
+
+// ----------------------------------------
+// Examining constant sqrt(2pi)
+//
+// Value bit pattern: 1.0100000011011001001100011111111101100010011100000101|10011 * 2^1
+// Value as double: 2.50662827463100068570156508940272033214569091796875
+// Relative error: -7.31205e-17
+// Correct digits: 16.14
+//
+// Inverse bit pattern: 1.1001100010000100010100110011110101000011011001010000|10001 * 2^-2
+// Inverse as double: 0.398942280401432702863218082711682654917240142822265625
+// Relative error: -6.24734e-17
+// Correct digits: 16.20
+// Can inverse be computed as double: no
+//
+// Conclusion: for sqrt(2pi), use inverse
+// ----------------------------------------
+
+
+
+static GOQuad qln10;
+
+static void
+print_bits (GOQuad const *qc)
+{
+       GOQuad qx = *qc;
+       GOQuad qhalf;
+       int e;
+       double s;
+       int b;
+       int N = 53 + 5;
+
+       qhalf.h = 0.5;
+       qhalf.l = 0;
+
+       if (qx.h < 0) {
+               g_printerr ("-");
+               qx.h = -qx.h;
+               qx.l = -qx.l;
+       }
+
+       (void)frexp (go_quad_value (&qx), &e);
+       s = ldexp (1.0, -e);
+       qx.h *= s;
+       qx.l *= s;
+
+       for (b = 0; b < N; b++) {
+               GOQuad d;
+               qx.h *= 2;
+               qx.l *= 2;
+               // For the last bit we compare against 1/2 to simulate rounding.
+               go_quad_sub (&d, &qx, (b == (N - 1) ? &qhalf : &go_quad_one));
+               if (go_quad_value (&d) >= 0) {
+                       qx = d;
+                       g_printerr ("1");
+               } else {
+                       g_printerr ("0");
+               }
+               if (b == 0)
+                       g_printerr (".");
+               if (b == 52)
+                       g_printerr ("|");
+       }
+
+       g_printerr (" * 2^%d", e - 1);
+}
+
+
+static double
+qmlogabs (GOQuad const *qx)
+{
+       GOQuad qy;
+       qy.h = fabs (qx->h); qy.l = fabs (qx->l);
+       go_quad_log (&qy, &qy);
+       go_quad_div (&qy, &qy, &qln10);
+       return -go_quad_value (&qy);
+}
+
+
+static void
+examine_constant (const char *descr, GOQuad const *qc)
+{
+       double dc, dic, dcic;
+       GOQuad qd, qic;
+       double dig_direct, dig_inverse;
+
+       g_printerr ("-----------------------------------------------------------------------------\n");
+       g_printerr ("Examining constant %s\n", descr);
+       g_printerr ("\n");
+
+       dc = go_quad_value (qc);
+       g_printerr ("Value bit pattern: ");
+       print_bits (qc);
+       g_printerr ("\n");
+       g_printerr ("Value as double: %.55g\n", dc);
+       go_quad_init (&qd, dc);
+       go_quad_sub (&qd, qc, &qd);
+       go_quad_div (&qd, &qd, qc);
+       g_printerr ("Relative error: %g\n", go_quad_value (&qd));
+       dig_direct = qmlogabs (&qd);
+       g_printerr ("Correct digits: %.2f\n", dig_direct);
+
+       g_printerr ("\n");
+
+       // For the inverse we assume that using double-double is good enough to
+       // answer our question.
+       go_quad_div (&qic, &go_quad_one, qc);
+       dic = go_quad_value (&qic);
+       g_printerr ("Inverse bit pattern: ");
+       print_bits (&qic);
+       g_printerr ("\n");
+       g_printerr ("Inverse as double: %.55g\n", dic);
+       go_quad_init (&qd, dic);
+       go_quad_sub (&qd, &qic, &qd);
+       go_quad_div (&qd, &qd, &qic);
+       g_printerr ("Relative error: %g\n", go_quad_value (&qd));
+       dig_inverse = qmlogabs (&qd);
+       g_printerr ("Correct digits: %.2f\n", dig_inverse);
+
+       dcic = (double)(1 / dc);
+       g_printerr ("Can inverse be computed as double: %s\n",
+                   (dcic == dic ? "yes" : "no"));
+
+       g_printerr ("\n");
+
+       g_printerr ("Conclusion: for %s, use %s\n", descr,
+                   (dig_direct >= dig_inverse ? "direct" : "inverse"));
+}
+
+int
+main (int argc, char **argv)
+{
+       void *state;
+       GOQuad qc;
+
+       g_printerr ("Determine for certain constants, c, whether its representation\n");
+       g_printerr ("as a double is more accurate than its inverse's.\n\n");
+
+       state = go_quad_start ();
+
+       go_quad_init (&qln10, 10);
+       go_quad_log (&qln10, &qln10);
+
+       examine_constant ("pi", &go_quad_pi);
+       examine_constant ("e", &go_quad_e);
+       examine_constant ("EulerGamma", &go_quad_euler);
+
+       examine_constant ("log(2)", &go_quad_ln2);
+       examine_constant ("log(10)", &qln10);
+       go_quad_div (&qc, &go_quad_ln2, &qln10);
+       examine_constant ("log10(2)", &qc);
+
+       examine_constant ("sqrt(2)", &go_quad_sqrt2);
+       go_quad_init (&qc, 3);
+       go_quad_sqrt (&qc, &qc);
+       examine_constant ("sqrt(3)", &qc);
+       go_quad_init (&qc, 5);
+       go_quad_sqrt (&qc, &qc);
+       examine_constant ("sqrt(5)", &qc);
+       go_quad_sqrt (&qc, &go_quad_pi);
+       examine_constant ("sqrt(pi)", &qc);
+       go_quad_sqrt (&qc, &go_quad_2pi);
+       examine_constant ("sqrt(2pi)", &qc);
+
+       go_quad_init (&qc, 180);
+       go_quad_div (&qc, &qc, &go_quad_pi);
+       examine_constant ("180/pi", &qc);
+
+       go_quad_end (state);
+
+       return 0;
+}