[goffice] Matrix: minor accuracy improvement to matrix inversion.

[Morten Welinder <mortenw src gnome org>]

commit a5894366ac41d3efe040e86a944f7a514275ae69
Author: Morten Welinder <terra gnome org>
Date:   Tue Jul 24 20:21:12 2012 -0400

    Matrix: minor accuracy improvement to matrix inversion.

 ChangeLog                    |   10 ++++++++++
 NEWS                         |    3 +++
 goffice/math/go-quad.c       |    8 ++++++++
 goffice/math/go-regression.c |   26 ++++++++++++++++++++++----
 4 files changed, 43 insertions(+), 4 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 1918943..9d36997 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,13 @@
+2012-07-24  Morten Welinder  <terra gnome org>
+
+	* goffice/math/go-quad.c (go_quad_add): Complain if not properly
+	initialized.
+	(go_quad_start): Keep track of initialization depth.
+
+	* goffice/math/go-regression.c (go_matrix_invert): Ensure go_quad
+	has been properly initialized.  Don't re-use partial results after
+	they have been rounded to "double".
+
 2012-07-22  Jean Brefort  <jean brefort normalesup org>
 
 	* goffice/Makefile.am: make introspection stuff build when no libgoffice
diff --git a/NEWS b/NEWS
index a58070e..1cde308 100644
--- a/NEWS
+++ b/NEWS
@@ -1,5 +1,8 @@
 goffice 0.9.6:
 
+Morten:
+	* Minute accuracy improvement to matrix inversion.
+
 --------------------------------------------------------------------------
 goffice 0.9.5:
 
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index ba0e048..b0d6ab5 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -98,6 +98,8 @@ SUFFIX(go_quad_functional) (void)
 #endif
 }
 
+static guint SUFFIX(go_quad_depth) = 0;
+
 static DOUBLE SUFFIX(CST);
 
 void *
@@ -105,6 +107,9 @@ SUFFIX(go_quad_start) (void)
 {
 	void *res = NULL;
 
+	if (SUFFIX(go_quad_depth)++ > 0)
+		return NULL;
+
 	if (!SUFFIX(go_quad_functional) ())
 		g_warning ("quad precision math may not be completely accurate.");
 
@@ -133,6 +138,7 @@ SUFFIX(go_quad_start) (void)
 void
 SUFFIX(go_quad_end) (void *state)
 {
+	SUFFIX(go_quad_depth)--;
 	if (!state)
 		return;
 
@@ -167,6 +173,8 @@ SUFFIX(go_quad_add) (QUAD *res, const QUAD *a, const QUAD *b)
 		: b->h - r + a->h + a->l + b->l;
 	res->h = r + s;
 	res->l = r - res->h + s;
+
+	g_return_if_fail (SUFFIX(go_quad_depth) > 0);
 }
 
 void
diff --git a/goffice/math/go-regression.c b/goffice/math/go-regression.c
index 6728a8b..4f1f998 100644
--- a/goffice/math/go-regression.c
+++ b/goffice/math/go-regression.c
@@ -145,6 +145,19 @@ general_linear_regressionl (long double *const *const xss, int xdim,
 	}							\
   } while (0)
 
+#define PRINT_QMATRIX(var,dim1,dim2)				\
+  do {								\
+	int _i, _j, _d1, _d2;					\
+	_d1 = (dim1);						\
+	_d2 = (dim2);						\
+	for (_j = 0; _j < _d2; _j++) {				\
+	  for (_i = 0; _i < _d1; _i++) {			\
+		  g_printerr (" %12.8" FORMAT_g, SUFFIX(go_quad_value) (&(var)[_i][_j])); \
+	  }							\
+	  g_printerr ("\n");					\
+	}							\
+  } while (0)
+
 #endif
 
 /*
@@ -585,6 +598,7 @@ SUFFIX(go_matrix_invert) (MATRIX A, int n)
 	QMATRIX R;
 	GORegressionResult err;
 	gboolean has_result;
+	void *state = SUFFIX(go_quad_start) ();
 
 	ALLOC_MATRIX2 (Q, n, n, QUAD);
 	ALLOC_MATRIX2 (R, n, n, QUAD);
@@ -594,25 +608,29 @@ SUFFIX(go_matrix_invert) (MATRIX A, int n)
 
 	if (has_result) {
 		int i, j, k;
+		QUAD *x = g_new (QUAD, n);
+
 		for (k = 0; k < n; k++) {
 			/* Solve R x = Q^T e_k */
 			for (i = n - 1; i >= 0; i--) {
 				QUAD d = Q[i][k];
 				for (j = i + 1; j < n; j++) {
 					QUAD p;
-					SUFFIX(go_quad_init) (&p, A[k][j]);
-					SUFFIX(go_quad_mul) (&p, &R[j][i], &p);
+					SUFFIX(go_quad_mul) (&p, &R[j][i], &x[j]);
 					SUFFIX(go_quad_sub) (&d, &d, &p);
 				}
-				SUFFIX(go_quad_div) (&d, &d, &R[i][i]);
-				A[k][i] = SUFFIX(go_quad_value) (&d);
+				SUFFIX(go_quad_div) (&x[i], &d, &R[i][i]);
+				A[k][i] = SUFFIX(go_quad_value) (&x[i]);
 			}
 		}
+		g_free (x);
 	}
 
 	FREE_MATRIX (Q, n, n);
 	FREE_MATRIX (R, n, n);
 
+	SUFFIX(go_quad_end) (state);
+
 	return has_result;
 }