GNOME.org
 

[goffice] regression: improve accuracy a bit more.

commit c5bd02345bfbc79f8279a685831d89ad3150a1fe
Author: Morten Welinder <terra gnome org>
Date:   Sun May 16 09:06:07 2010 -0400

    regression: improve accuracy a bit more.

 ChangeLog                    |    5 +
 configure.in                 |    2 +
 goffice/math/Makefile.am     |   22 +++--
 goffice/math/go-quad.c       |  203 ++++++++++++++++++++++++++++++++++++++++++
 goffice/math/go-quad.h       |   52 +++++++++++
 goffice/math/go-regression.c |  127 ++++++++++++++++++--------
 goffice/math/goffice-math.h  |    1 +
 7 files changed, 364 insertions(+), 48 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 8687b68..64c0dda 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2010-05-16  Morten Welinder  <terra gnome org>
+
+	* goffice/math/go-quad.c: quad precision routines.
+	* goffice/math/go-regression.c (refine): Use quad precision.
+
 2010-05-13  Morten Welinder  <terra gnome org>
 
 	* goffice/component/go-component-factory.c (go_component_type_service_read_xml):
diff --git a/configure.in b/configure.in
index 5e68813..412e93c 100644
--- a/configure.in
+++ b/configure.in
@@ -394,6 +394,8 @@ AC_CHECK_HEADERS(ieeefp.h ieee754.h)
 dnl Check for implementations of passwd/group file entry functions
 AC_CHECK_HEADERS(pwd.h grp.h)
 
+AC_CHECK_HEADERS(fpu_control.h)
+
 dnl Check for some functions
 AC_CHECK_FUNCS(random drand48 finite memmove mkdtemp uname times sysconf)
 
diff --git a/goffice/math/Makefile.am b/goffice/math/Makefile.am
index 70a8bd1..4ca683e 100644
--- a/goffice/math/Makefile.am
+++ b/goffice/math/Makefile.am
@@ -2,12 +2,13 @@ noinst_LTLIBRARIES = libgoffice-math.la
 
 libgoffice_math_la_SOURCES =	\
 	go-math.c		\
-	go-rangefunc.c	\
-	go-regression.c	\
-	go-cspline.c	\
-	go-complex.c	\
+	go-rangefunc.c		\
+	go-regression.c		\
+	go-cspline.c		\
+	go-complex.c		\
 	go-fft.c		\
-	go-matrix3x3.c	\
+	go-matrix3x3.c		\
+	go-quad.c		\
 	go-R.c			\
 	go-distribution.c
 
@@ -15,12 +16,13 @@ libgoffice_math_ladir = $(goffice_include_dir)/math
 libgoffice_math_la_HEADERS = 	\
 	goffice-math.h		\
 	go-math.h		\
-	go-rangefunc.h	\
-	go-regression.h	\
-	go-cspline.h	\
-	go-complex.h	\
+	go-rangefunc.h		\
+	go-regression.h		\
+	go-cspline.h		\
+	go-complex.h		\
 	go-fft.h		\
-	go-matrix3x3.h	\
+	go-matrix3x3.h		\
+	go-quad.h		\
 	go-R.h			\
 	go-distribution.h
 
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
new file mode 100644
index 0000000..eeb0aef
--- /dev/null
+++ b/goffice/math/go-quad.c
@@ -0,0 +1,203 @@
+/*
+ * go-quad.c:  Extended precision routines.
+ *
+ * Authors:
+ *   Morten Welinder <terra gnome org>
+ *
+ * This follows "A Floating-Point Technoque for Extending the Available
+ * Precision" by T. J. Dekker in _Numerische Mathematik_ 18.
+ * Springer Verlag 1971.
+ *
+ * Note: for this to work, the processor must evaluate with the right
+ * precision.  For ix86 that means trouble as the default is to evaluate
+ * with long-double precision internally.  We solve this by setting the
+ * relevant x87 control flag.
+ *
+ * Note: the compiler should not rearrange expressions.
+ */
+
+#include <goffice/goffice-config.h>
+#include "go-quad.h"
+#include <math.h>
+
+/* Normalize cpu id.  */
+#if !defined(i386) && (defined(__i386__) || defined(__i386))
+#define i386 1
+#endif
+
+#if defined(i386) && defined(HAVE_FPU_CONTROL_H)
+#include <fpu_control.h>
+#endif
+
+#ifndef DOUBLE
+
+#define QUAD SUFFIX(GOQuad)
+
+#define DOUBLE double
+#define SUFFIX(_n) _n
+#define DOUBLE_MANT_DIG DBL_MANT_DIG
+
+#ifdef GOFFICE_WITH_LONG_DOUBLE
+#include "go-quad.c"
+#undef DOUBLE
+#undef SUFFIX
+#undef DOUBLE_MANT_DIG
+#define DOUBLE long double
+#define SUFFIX(_n) _n ## l
+#define DOUBLE_MANT_DIG LDBL_MANT_DIG
+#endif
+
+#endif
+
+gboolean
+SUFFIX(go_quad_functional) (void)
+{
+	if (FLT_RADIX != 2)
+		return FALSE;
+
+#ifdef i386
+	if (sizeof (DOUBLE) != sizeof (double))
+		return TRUE;
+
+#ifdef HAVE_FPU_CONTROL_H
+	return TRUE;
+#else
+	return FALSE;
+#endif
+#else
+	return TRUE;
+#endif
+}
+
+static DOUBLE SUFFIX(CST);
+
+void *
+SUFFIX(go_quad_start) (void)
+{
+	void *res = NULL;
+
+	if (!go_quad_functional ())
+		g_warning ("quad precision math may not be completely accurate.");
+
+#ifdef i386
+	if (sizeof (DOUBLE) == sizeof (double)) {
+#ifdef HAVE_FPU_CONTROL_H
+		fpu_control_t state, newstate;
+		fpu_control_t mask =
+			_FPU_EXTENDED | _FPU_DOUBLE | _FPU_SINGLE;
+		_FPU_GETCW (state);
+		newstate = (state & ~mask) | _FPU_DOUBLE;
+		_FPU_SETCW (newstate);
+
+		res = g_memdup (&state, sizeof (state));
+#else
+		/* Hope for the best.  */
+#endif
+	}
+#endif
+
+	SUFFIX(CST) = 1 + SUFFIX(ldexp) (1.0, (DOUBLE_MANT_DIG + 1) / 2);
+
+	return res;
+}
+
+void
+SUFFIX(go_quad_end) (void *state)
+{
+	if (!state)
+		return;
+
+#ifdef i386
+#ifdef HAVE_FPU_CONTROL_H
+	_FPU_SETCW (*(fpu_control_t*)state);
+#endif
+#endif
+
+	g_free (state);
+}
+
+void
+SUFFIX(go_quad_init) (QUAD *res, DOUBLE h, DOUBLE l)
+{
+	res->h = h;
+	res->l = l;
+}
+
+DOUBLE
+SUFFIX(go_quad_value) (const QUAD *a)
+{
+	return a->h + a->l;
+}
+
+void
+SUFFIX(go_quad_add) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+	DOUBLE r = a->h + b->h;
+	DOUBLE s = SUFFIX(fabs) (a->h) > SUFFIX(fabs) (b->h)
+		? a->h - r + b->h + b->l + a->l
+		: b->h - r + a->h + a->l + b->l;
+	res->h = r + s;
+	res->l = r - res->h + s;
+}
+
+void
+SUFFIX(go_quad_sub) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+	DOUBLE r = a->h - b->h;
+	DOUBLE s = SUFFIX(fabs) (a->h) > SUFFIX(fabs) (b->h)
+		? +a->h - r - b->h - b->l + a->l
+		: -b->h - r + a->h + a->l - b->l;
+	res->h = r + s;
+	res->l = r - res->h + s;
+}
+
+void
+SUFFIX(go_quad_mul12) (QUAD *res, DOUBLE x, DOUBLE y)
+{
+	DOUBLE p1 = x * SUFFIX(CST);
+	DOUBLE hx = x - p1 + p1;
+	DOUBLE tx = x - hx;
+
+	DOUBLE p2 = y * SUFFIX(CST);
+	DOUBLE hy = y - p2 + p2;
+	DOUBLE ty = y - hy;
+
+	DOUBLE p = hx * hy;
+	DOUBLE q = hx * ty + tx * hy;
+	res->h = p + q;
+	res->l = p - res->h + q + tx * ty;
+}
+
+void
+SUFFIX(go_quad_mul) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+	QUAD c;
+	SUFFIX(go_quad_mul12) (&c, a->h, b->h);
+	c.l = a->h * b->l + a->l * b->h + c.l;
+	res->h = c.h + c.l;
+	res->l = c.h - res->h + c.l;
+}
+
+void
+SUFFIX(go_quad_div) (QUAD *res, const QUAD *a, const QUAD *b)
+{
+	QUAD c, u;
+	c.h = a->h / b->h;
+	SUFFIX(go_quad_mul12) (&u, c.h, b->h);
+	c.l = (a->h - u.h - u.l + a->l - c.h * b->l) / b->h;
+	res->h = c.h + c.l;
+	res->l = c.h - res->h + c.l;
+}
+
+int
+SUFFIX(go_quad_sgn) (const QUAD *a)
+{
+	DOUBLE d = (a->h == 0 ? a->l : a->h);
+
+	if (d == 0)
+		return 0;
+	else if (d > 0)
+		return +1;
+	else
+		return -1;
+}
diff --git a/goffice/math/go-quad.h b/goffice/math/go-quad.h
new file mode 100644
index 0000000..53acf31
--- /dev/null
+++ b/goffice/math/go-quad.h
@@ -0,0 +1,52 @@
+#ifndef __GO_QUAD_H
+#define __GO_QUAD_H
+
+#include <glib.h>
+
+G_BEGIN_DECLS
+
+typedef struct {
+	double h;
+	double l;
+} GOQuad;
+
+gboolean go_quad_functional (void);
+void *go_quad_start (void);
+void go_quad_end (void *state);
+
+void go_quad_init (GOQuad *res, double h, double l);
+
+double go_quad_value (const GOQuad *a);
+void go_quad_add (GOQuad *res, const GOQuad *a, const GOQuad *b);
+void go_quad_sub (GOQuad *res, const GOQuad *a, const GOQuad *b);
+void go_quad_mul (GOQuad *res, const GOQuad *a, const GOQuad *b);
+void go_quad_div (GOQuad *res, const GOQuad *a, const GOQuad *b);
+int  go_quad_sgn (const GOQuad *a);
+
+void go_quad_mul12 (GOQuad *res, double x, double y);
+
+#ifdef GOFFICE_WITH_LONG_DOUBLE
+typedef struct {
+	long double h;
+	long double l;
+} GOQuadl;
+
+gboolean go_quad_functionall (void);
+void *go_quad_startl (void);
+void go_quad_endl (void *state);
+
+void go_quad_initl (GOQuadl *res, long double h, long double l);
+
+long double go_quad_valuel (const GOQuadl *a);
+void go_quad_addl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+void go_quad_subl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+void go_quad_mull (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+void go_quad_divl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
+int  go_quad_sgnl (const GOQuadl *a);
+
+void go_quad_mul12l (GOQuadl *res, long double x, long double y);
+#endif
+
+G_END_DECLS
+
+#endif
diff --git a/goffice/math/go-regression.c b/goffice/math/go-regression.c
index 006dd9f..45464f8 100644
--- a/goffice/math/go-regression.c
+++ b/goffice/math/go-regression.c
@@ -17,20 +17,21 @@
 #include <stdlib.h>
 #include <stdio.h>
 
+/* ------------------------------------------------------------------------- */
+/* Self-inclusion magic.  */
+
 #ifndef DOUBLE
 
 #define DEFINE_COMMON
 #define DOUBLE double
 #define DOUBLE_MANT_DIG DBL_MANT_DIG
 #define SUFFIX(_n) _n
+#define FORMAT_f "f"
 #define FORMAT_g "g"
 
-#define MATRIX DOUBLE **
-#define CONSTMATRIX DOUBLE *const *const
-
-#if 0
 #ifdef GOFFICE_WITH_LONG_DOUBLE
 
+#if 0
 /* We define this to cause certain "double" functions to be implemented in
    terms of their "long double" counterparts.  */
 #define BOUNCE
@@ -49,6 +50,7 @@ general_linear_regressionl (long double *const *const xss, int xdim,
 #undef DOUBLE
 #undef DOUBLE_MANT_DIG
 #undef SUFFIX
+#undef FORMAT_f
 #undef FORMAT_g
 #ifdef HAVE_SUNMATH_H
 #include <sunmath.h>
@@ -56,13 +58,20 @@ general_linear_regressionl (long double *const *const xss, int xdim,
 #define DOUBLE long double
 #define DOUBLE_MANT_DIG LDBL_MANT_DIG
 #define SUFFIX(_n) _n ## l
+#define FORMAT_f "Lf"
 #define FORMAT_g "Lg"
 #endif
 
 #endif
 
+/* ------------------------------------------------------------------------- */
+
 #ifdef DEFINE_COMMON
 
+#define MATRIX DOUBLE **
+#define CONSTMATRIX DOUBLE *const *const
+#define QUAD SUFFIX(GOQuad)
+
 #undef DEBUG_NEAR_SINGULAR
 #undef DEBUG_REFINEMENT
 
@@ -112,25 +121,25 @@ general_linear_regressionl (long double *const *const xss, int xdim,
 	int _i, _j, _d1, _d2;					\
 	_d1 = (dim1);						\
 	_d2 = (dim2);						\
-	for (_i = 0; _i < _d1; _i++)				\
-	  {							\
-	    for (_j = 0; _j < _d2; _j++)			\
-	      g_printerr (" %19.10" FORMAT_g, (var)[_i][_j]);	\
-	    g_printerr ("\n");					\
+	for (_j = 0; _j < _d2; _j++) {				\
+	  for (_i = 0; _i < _d1; _i++) {			\
+	    g_printerr (" %12.8" FORMAT_g, (var)[_i][_j]);	\
 	  }							\
+	  g_printerr ("\n");					\
+	}							\
   } while (0)
 
 #endif
 
 /*
- *       ---> j
+ *       ---> i
  *
  *  |    ********
  *  |    ********
  *  |    ********        A[i][j]
  *  v    ********
  *       ********
- *  i    ********
+ *  j    ********
  *       ********
  *       ********
  *
@@ -207,6 +216,10 @@ SUFFIX(QR) (CONSTMATRIX A, MATRIX Q, MATRIX R, int m, int n)
 
 		(void)SUFFIX(go_range_sumsq) (Q[k], n, &L);
 		L = SUFFIX(sqrt) (L);
+#if 0
+		PRINT_MATRIX (Q, m, n);
+		g_printerr ("L[%d] = %20.15" FORMAT_g "\n", k, L);
+#endif
 		if (L == 0)
 			return GO_REG_singular;
 
@@ -233,20 +246,28 @@ SUFFIX(QR) (CONSTMATRIX A, MATRIX Q, MATRIX R, int m, int n)
 
 static void
 SUFFIX(calc_residual) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
-		       const DOUBLE *res, DOUBLE *residual, DOUBLE *N2)
+		       const DOUBLE *res, QUAD *residual, QUAD *N2)
 {
 	int i, j;
 
-	*N2 = 0;
+	SUFFIX(go_quad_init) (N2, 0, 0);
 
 	for (i = 0; i < n; i++) {
-		DOUBLE d = 0;
-		for (j = 0; j < dim; j++)
-			d += A[j][i] * res[j];
-		d = b[i] - d;
-		(*N2) += d * d;
+		QUAD d;
+
+		SUFFIX(go_quad_init) (&d, b[i], 0);
+
+		for (j = 0; j < dim; j++) {
+			QUAD e;
+			SUFFIX(go_quad_mul12) (&e, A[j][i], res[j]);
+			SUFFIX(go_quad_sub) (&d, &d, &e);
+		}
+
 		if (residual)
 			residual[i] = d;
+
+		SUFFIX(go_quad_mul) (&d, &d, &d);
+		SUFFIX(go_quad_add) (N2, N2, &d);
 	}
 }
 
@@ -254,50 +275,79 @@ static void
 SUFFIX(refine) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
 		CONSTMATRIX Q, CONSTMATRIX R, DOUBLE *result)
 {
-	DOUBLE *residual = g_new (DOUBLE, n);
+	QUAD *residual = g_new (QUAD, n);
+	QUAD *newresidual = g_new (QUAD, n);
 	DOUBLE *delta = g_new (DOUBLE, dim);
 	DOUBLE *newresult = g_new (DOUBLE, dim);
 	int pass;
-	DOUBLE best_N2;
+	QUAD best_N2;
+	void *state;
+
+	state = SUFFIX(go_quad_start) ();
 
 	SUFFIX(calc_residual) (A, b, dim, n, result, residual, &best_N2);
 #ifdef DEBUG_REFINEMENT
-	g_printerr ("-: Residual norm = %20.15" FORMAT_g "\n", best_N2);
+	g_printerr ("-: Residual norm = %20.15" FORMAT_g "\n",
+		    SUFFIX(go_quad_value) (&best_N2));
 #endif
 
-	for (pass = 1; pass < 5; pass++) {
+	for (pass = 1; pass < 10; pass++) {
 		int i, j;
-		DOUBLE N2;
+		QUAD dres, N2;
 
 		/* newresult = R^-1 Q^T residual */
 		for (i = dim - 1; i >= 0; i--) {
-			DOUBLE acc = SUFFIX(dot_product) (Q[i], residual, n);
-			for (j = i + 1; j < dim; j++)
-				acc -= R[j][i] * delta[j];
-			delta[i] = acc / R[i][i];
+			QUAD acc, d, Rii;
+
+			SUFFIX(go_quad_init) (&acc, 0, 0);
+
+			for (j = 0; j < n; j++) {
+				QUAD q, qr;
+				SUFFIX(go_quad_init) (&q, Q[i][j], 0);
+				SUFFIX(go_quad_mul) (&qr, &q, residual + j);
+				SUFFIX(go_quad_add) (&acc, &acc, &qr);
+			}
+			for (j = i + 1; j < dim; j++) {
+				QUAD Rd;
+				SUFFIX(go_quad_mul12) (&Rd, R[j][i], delta[j]);
+				SUFFIX(go_quad_sub) (&acc, &acc, &Rd);
+			}
+
+			SUFFIX(go_quad_init) (&Rii, R[i][i], 0);
+			SUFFIX(go_quad_div) (&d, &acc, &Rii);
+
+			delta[i] = SUFFIX(go_quad_value) (&d);
+			newresult[i] = delta[i] + result[i];
 #ifdef DEBUG_REFINEMENT
-			g_printerr ("d[%2d] = %20" FORMAT_g
-				    "  %20" FORMAT_g "\n",
-				    i, delta[i], result[i]);
+			g_printerr ("d[%2d] = %20.15" FORMAT_f
+				    "  %20.15" FORMAT_f
+				    "  %20.15" FORMAT_f "\n",
+				    i, delta[i], result[i], newresult[i]);
 #endif
 		}
 
-		for (i = 0; i < dim; i++)
-			newresult[i] = delta[i] + result[i];
-
-		SUFFIX(calc_residual) (A, b, dim, n, newresult, residual, &N2);
+		SUFFIX(calc_residual) (A, b, dim, n, newresult, newresidual, &N2);
+		SUFFIX (go_quad_sub) (&dres, &N2, &best_N2);
 
 #ifdef DEBUG_REFINEMENT
-		g_printerr ("%d: Residual norm = %20.15" FORMAT_g "\n", pass, N2);
+		g_printerr ("%d: Residual norm = %20.15" FORMAT_g
+			    " (%.5" FORMAT_g ")\n",
+			    pass,
+			    SUFFIX(go_quad_value) (&N2),
+			    SUFFIX(go_quad_value) (&dres));
 #endif
-		if (N2 >= best_N2)
+		if (SUFFIX(go_quad_sgn) (&dres) >= 0)
 			break;
 
 		best_N2 = N2;
 		COPY_VECTOR (result, newresult, dim);
+		COPY_VECTOR (residual, newresidual, n);
 	}
 
+	SUFFIX(go_quad_end) (state);
+
 	g_free (residual);
+	g_free (newresidual);
 	g_free (newresult);
 	g_free (delta);
 }
@@ -661,6 +711,7 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
 		DOUBLE *inv = g_new (DOUBLE, xdim);
 		int err;
 		int k;
+		QUAD N2;
 
 		/* This should not fail since n >= 1.  */
 		err = SUFFIX(go_range_average) (ys, n, &stat_->ybar);
@@ -680,8 +731,8 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
 			g_assert (err == 0);
 		}
 
-		SUFFIX(calc_residual) (xss, ys, xdim, n, result, NULL,
-				       &stat_->ss_resid);
+		SUFFIX(calc_residual) (xss, ys, xdim, n, result, NULL, &N2);
+		stat_->ss_resid = SUFFIX(go_quad_value) (&N2);
 
 		stat_->sqr_r = (stat_->ss_total == 0)
 			? 1
diff --git a/goffice/math/goffice-math.h b/goffice/math/goffice-math.h
index 7762d3c..c307224 100644
--- a/goffice/math/goffice-math.h
+++ b/goffice/math/goffice-math.h
@@ -9,6 +9,7 @@
 #include <goffice/math/go-matrix3x3.h>
 #include <goffice/math/go-rangefunc.h>
 #include <goffice/math/go-regression.h>
+#include <goffice/math/go-quad.h>
 #include <goffice/math/go-R.h>
 
 #endif
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]