GNOME.org
 

[goffice] Rangefuncs: improve precision

commit 8ea33770839d1132647b8d3d056706c06b140ccb
Author: Morten Welinder <terra gnome org>
Date:   Wed Nov 23 15:23:23 2011 -0500

    Rangefuncs: improve precision

 ChangeLog                     |    8 ++
 NEWS                          |    1 +
 goffice/math/Makefile.am      |    2 +
 goffice/math/go-accumulator.c |  134 ++++++++++++++++++++++++++++++++++
 goffice/math/go-accumulator.h |   38 ++++++++++
 goffice/math/go-quad.c        |    2 +-
 goffice/math/go-rangefunc.c   |  160 ++++++++++++++++++++++-------------------
 goffice/math/go-rangefunc.h   |    2 +
 goffice/math/goffice-math.h   |    5 +-
 9 files changed, 274 insertions(+), 78 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index fa2d115..da2ec57 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+2011-11-23  Morten Welinder  <terra gnome org>
+
+	* goffice/math/go-accumulator.c: New file.
+
+	* goffice/math/go-rangefunc.c (go_range_devsq, go_range_sumsq)
+	(go_range_sum, go_range_average): Base on GOAccumulator.
+	(go_range_constant): New function.
+
 2011-10-27  Andreas J. Guelzow <aguelzow pyrshep ca>
 
 	* goffice/gtk/goffice-gtk.c (recent_func_log_func): new
diff --git a/NEWS b/NEWS
index 586ffbd..8c9b22f 100644
--- a/NEWS
+++ b/NEWS
@@ -48,6 +48,7 @@ Morten:
 	* Recognize scientific formats with longer exponents.
 	* Fix problem with sticky format colouring.
 	* Improve certain range functions in corner cases.  [#660564]
+	* Make go_range_sum perfect, modulo overflow.  [#664656]
 	* Add resource manager for embedded files.
 	* Plug leaks.
 
diff --git a/goffice/math/Makefile.am b/goffice/math/Makefile.am
index 4ca683e..94f2c73 100644
--- a/goffice/math/Makefile.am
+++ b/goffice/math/Makefile.am
@@ -1,6 +1,7 @@
 noinst_LTLIBRARIES = libgoffice-math.la
 
 libgoffice_math_la_SOURCES =	\
+	go-accumulator.c	\
 	go-math.c		\
 	go-rangefunc.c		\
 	go-regression.c		\
@@ -15,6 +16,7 @@ libgoffice_math_la_SOURCES =	\
 libgoffice_math_ladir = $(goffice_include_dir)/math
 libgoffice_math_la_HEADERS = 	\
 	goffice-math.h		\
+	go-accumulator.h	\
 	go-math.h		\
 	go-rangefunc.h		\
 	go-regression.h		\
diff --git a/goffice/math/go-accumulator.c b/goffice/math/go-accumulator.c
new file mode 100644
index 0000000..a9d44aa
--- /dev/null
+++ b/goffice/math/go-accumulator.c
@@ -0,0 +1,134 @@
+/*
+ * go-accumulator.c: high-precision sums
+ *
+ * Authors:
+ *   Morten Welinder <terra gnome org>
+ *
+ * See http://code.activestate.com/recipes/393090/
+ *
+ * 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-accumulator.h"
+#include "go-quad.h"
+#include <math.h>
+#include <float.h>
+
+#ifndef DOUBLE
+
+#define DOUBLE double
+#define SUFFIX(_n) _n
+#define DOUBLE_MANT_DIG DBL_MANT_DIG
+
+struct GOAccumulator_ {
+	GArray *partials;
+};
+#define ACC SUFFIX(GOAccumulator)
+
+#ifdef GOFFICE_WITH_LONG_DOUBLE
+#include "go-accumulator.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
+
+struct GOAccumulatorl_ {
+	GArray *partials;
+};
+
+#endif
+#endif
+
+
+
+
+gboolean
+SUFFIX(go_accumulator_functional) (void)
+{
+	return SUFFIX(go_quad_functional) ();
+}
+
+void *
+SUFFIX(go_accumulator_start) (void)
+{
+	return SUFFIX(go_quad_start) ();
+}
+
+void
+SUFFIX(go_accumulator_end) (void *state)
+{
+	SUFFIX(go_quad_end) (state);
+}
+
+ACC *
+SUFFIX(go_accumulator_new) (void)
+{
+	ACC *res = g_new (ACC, 1);
+	res->partials = g_array_sized_new (FALSE, FALSE, sizeof (DOUBLE), 10);
+	return res;
+}
+
+void
+SUFFIX(go_accumulator_free) (ACC *acc)
+{
+	g_return_if_fail (acc != NULL);
+
+	g_array_free (acc->partials, TRUE);
+	g_free (acc);
+}
+
+void
+SUFFIX(go_accumulator_clear) (ACC *acc)
+{
+	g_return_if_fail (acc != NULL);
+	g_array_set_size (acc->partials, 0);
+}
+
+void
+SUFFIX(go_accumulator_add) (ACC *acc, DOUBLE x)
+{
+	unsigned ui = 0, uj;
+
+	g_return_if_fail (acc != NULL);
+
+	for (uj = 0; uj < acc->partials->len; uj++) {
+		DOUBLE y = g_array_index (acc->partials, DOUBLE, uj);
+		DOUBLE hi, lo;
+		if (SUFFIX(fabs)(x) < SUFFIX(fabs)(y)) {
+			DOUBLE t = x;
+			x = y;
+			y = t;
+		}
+		hi = x + y;
+		lo = y - (hi - x);
+		if (lo != 0) {
+			g_array_index (acc->partials, DOUBLE, ui) = lo;
+			ui++;
+		}
+		x = hi;
+	}
+	g_array_set_size (acc->partials, ui + 1);
+	g_array_index (acc->partials, DOUBLE, ui) = x;
+}
+
+DOUBLE
+SUFFIX(go_accumulator_value) (ACC *acc)
+{
+	unsigned ui;
+	DOUBLE sum = 0;
+
+	g_return_val_if_fail (acc != NULL, 0);
+
+	for (ui = 0; ui < acc->partials->len; ui++)
+		sum += g_array_index (acc->partials, DOUBLE, ui);
+
+	return sum;
+}
diff --git a/goffice/math/go-accumulator.h b/goffice/math/go-accumulator.h
new file mode 100644
index 0000000..1b25d96
--- /dev/null
+++ b/goffice/math/go-accumulator.h
@@ -0,0 +1,38 @@
+#ifndef __GO_ACCUMULATOR_H
+#define __GO_ACCUMULATOR_H
+
+#include <glib.h>
+
+G_BEGIN_DECLS
+
+typedef struct GOAccumulator_ GOAccumulator;
+
+gboolean go_accumulator_functional (void);
+void *go_accumulator_start (void);
+void go_accumulator_end (void *state);
+
+GOAccumulator *go_accumulator_new (void);
+void go_accumulator_free (GOAccumulator *acc);
+void go_accumulator_clear (GOAccumulator *acc);
+void go_accumulator_add (GOAccumulator *acc, double x);
+double go_accumulator_value (GOAccumulator *acc);
+
+
+#ifdef GOFFICE_WITH_LONG_DOUBLE
+
+typedef struct GOAccumulatorl_ GOAccumulatorl;
+
+gboolean go_accumulator_functionall (void);
+void *go_accumulator_startl (void);
+void go_accumulator_endl (void *state);
+
+GOAccumulatorl *go_accumulator_newl (void);
+void go_accumulator_freel (GOAccumulatorl *acc);
+void go_accumulator_clearl (GOAccumulatorl *acc);
+void go_accumulator_addl (GOAccumulatorl *acc, long double x);
+long double go_accumulator_valuel (GOAccumulatorl *acc);
+#endif
+
+G_END_DECLS
+
+#endif
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index 8d76934..ba87492 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -4,7 +4,7 @@
  * Authors:
  *   Morten Welinder <terra gnome org>
  *
- * This follows "A Floating-Point Technoque for Extending the Available
+ * This follows "A Floating-Point Technique for Extending the Available
  * Precision" by T. J. Dekker in _Numerische Mathematik_ 18.
  * Springer Verlag 1971.
  *
diff --git a/goffice/math/go-rangefunc.c b/goffice/math/go-rangefunc.c
index 1faa6f2..70e2eae 100644
--- a/goffice/math/go-rangefunc.c
+++ b/goffice/math/go-rangefunc.c
@@ -8,6 +8,8 @@
 
 #include <goffice/goffice-config.h>
 #include "go-rangefunc.h"
+#include "go-accumulator.h"
+#include "go-quad.h"
 
 #include <math.h>
 #include <stdlib.h>
@@ -42,54 +44,27 @@
 
 /* ------------------------------------------------------------------------- */
 
-static int
-SUFFIX(identical_helper) (DOUBLE const *xs, int n)
+static SUFFIX(GOAccumulator)*
+SUFFIX(sum_helper) (DOUBLE const *xs, int n)
 {
-	DOUBLE x;
-	int i;
-
-	if (n <= 1)
-		return n;
-
-	x = xs[0];
-	for (i = 1; i < n; i++)
-		if (x != xs[i])
-			break;
-
-	return i;
-}
-
-static LDOUBLE
-SUFFIX(sum_helper) (DOUBLE const *xs, int n, int *all_id)
-{
-	LDOUBLE sum;
-	int i;
-
-	/*
-	 * We treat the an initial constant part special in order to catch
-	 * the important special case of completely constant.  We do not
-	 * want to accumulate rounding errors even though the extra 11?
-	 * bits in long double will shield us from n<2^11, give or take.
-	 */
-
-	i = SUFFIX(identical_helper) (xs, n);
-	*all_id = (i == n);
-	sum = i ? i * (LDOUBLE)(xs[0]) : 0;
-
-	for (; i < n; i++)
-		sum += xs[i];
-
-	return sum;
+	SUFFIX(GOAccumulator) *acc = SUFFIX(go_accumulator_new) ();
+	while (n > 0) {
+		n--;
+		SUFFIX(go_accumulator_add) (acc, xs[n]);
+	}
+	return acc;
 }
 
-/* ------------------------------------------------------------------------- */
 
 /* Arithmetic sum.  */
 int
 SUFFIX(go_range_sum) (DOUBLE const *xs, int n, DOUBLE *res)
 {
-	int all_id;
-	*res = SUFFIX(sum_helper) (xs, n, &all_id);
+	void *state = SUFFIX(go_accumulator_start) ();
+	SUFFIX(GOAccumulator) *acc = SUFFIX(sum_helper) (xs, n);
+	*res = SUFFIX(go_accumulator_value) (acc);
+	SUFFIX(go_accumulator_free) (acc);
+	SUFFIX(go_accumulator_end) (state);
 	return 0;
 }
 
@@ -97,20 +72,19 @@ SUFFIX(go_range_sum) (DOUBLE const *xs, int n, DOUBLE *res)
 int
 SUFFIX(go_range_sumsq) (DOUBLE const *xs, int n, DOUBLE *res)
 {
-	/* http://bugzilla.gnome.org/show_bug.cgi?id=131588 */
-#ifdef HAVE_LONG_DOUBLE
-	long double x, sum = 0;
-#else
-	DOUBLE x, sum = 0;
-#endif
-	int i;
-
-	for (i = 0; i < n; i++) {
-		x = xs[i];
-		sum += x * x;
+	void *state = SUFFIX(go_accumulator_start) ();
+	SUFFIX(GOAccumulator) *acc = SUFFIX(go_accumulator_new) ();
+	while (n > 0) {
+		SUFFIX(GOQuad) q;
+		n--;
+		SUFFIX(go_quad_init) (&q, xs[n]);
+		SUFFIX(go_quad_mul) (&q, &q, &q);
+		SUFFIX(go_accumulator_add) (acc, q.h);
+		SUFFIX(go_accumulator_add) (acc, q.l);
 	}
-
-	*res = sum;
+	*res = SUFFIX(go_accumulator_value) (acc);
+	SUFFIX(go_accumulator_free) (acc);
+	SUFFIX(go_accumulator_end) (state);
 	return 0;
 }
 
@@ -118,15 +92,16 @@ SUFFIX(go_range_sumsq) (DOUBLE const *xs, int n, DOUBLE *res)
 int
 SUFFIX(go_range_average) (DOUBLE const *xs, int n, DOUBLE *res)
 {
-	int all_id;
-
 	if (n <= 0)
 		return 1;
 
-	*res = SUFFIX(sum_helper) (xs, n, &all_id) / n;
-	if (all_id)
+	if (SUFFIX(go_range_constant) (xs, n)) {
 		*res = xs[0];
+		return 0;
+	}
 
+	SUFFIX(go_range_sum) (xs, n, res);
+	*res /= n;
 	return 0;
 }
 
@@ -172,9 +147,11 @@ SUFFIX(go_range_maxabs) (DOUBLE const *xs, int n, DOUBLE *res)
 		DOUBLE max = SUFFIX(fabs) (xs[0]);
 		int i;
 
-		for (i = 1; i < n; i++)
-			if (SUFFIX(fabs) (xs[i]) > max)
-				max = SUFFIX(fabs) (xs[i]);
+		for (i = 1; i < n; i++) {
+			DOUBLE xabs = SUFFIX(fabs) (xs[i]);
+			if (xabs > max)
+				max = xabs;
+		}
 		*res = max;
 		return 0;
 	} else
@@ -185,23 +162,46 @@ SUFFIX(go_range_maxabs) (DOUBLE const *xs, int n, DOUBLE *res)
 int
 SUFFIX(go_range_devsq) (DOUBLE const *xs, int n, DOUBLE *res)
 {
-	LDOUBLE q = 0;
-
-	if (n > 0) {
-		int i, all_id;
-		LDOUBLE m = SUFFIX(sum_helper) (xs, n, &all_id) / n;
-
-		if (all_id) {
-			*res = 0;
-			return 0;
-		}
-
-		for (i = 0; i < n; i++) {
-			LDOUBLE dx = xs[i] - m;
-			q += dx * dx;
+	if (SUFFIX(go_range_constant) (xs, n))
+		*res = 0;
+	else {
+		void *state;
+		SUFFIX(GOAccumulator) *acc;
+		DOUBLE sum, sumh, suml;
+		SUFFIX(GOQuad) qavg, qtmp, qn;
+
+		state = SUFFIX(go_accumulator_start) ();
+		acc = SUFFIX(sum_helper) (xs, n);
+		sumh = SUFFIX(go_accumulator_value) (acc);
+		SUFFIX(go_accumulator_add) (acc, -sumh);
+		suml = SUFFIX(go_accumulator_value) (acc);
+
+		SUFFIX(go_quad_init) (&qavg, sumh);
+		SUFFIX(go_quad_init) (&qtmp, suml);
+		SUFFIX(go_quad_add) (&qavg, &qavg, &qtmp);
+		SUFFIX(go_range_sum) (xs, n, &sum);
+		SUFFIX(go_quad_init) (&qn, n);
+		SUFFIX(go_quad_div) (&qavg, &qavg, &qn);
+		/*
+		 * The just-generated qavg isn't exact, but should be
+		 * good to near-GOQuad precision.  We hope that's good
+		 * enough.
+		 */
+
+		SUFFIX(go_accumulator_clear) (acc);
+		while (n > 0) {
+			SUFFIX(GOQuad) q;
+			n--;
+			SUFFIX(go_quad_init) (&q, xs[n]);
+			SUFFIX(go_quad_sub) (&q, &q, &qavg);
+			SUFFIX(go_quad_mul) (&q, &q, &q);
+			SUFFIX(go_accumulator_add) (acc, q.h);
+			SUFFIX(go_accumulator_add) (acc, q.l);
 		}
+		*res = SUFFIX(go_accumulator_value) (acc);
+		SUFFIX(go_accumulator_free) (acc);
+		SUFFIX(go_accumulator_end) (state);
 	}
-	*res = q;
 	return 0;
 }
 
@@ -285,6 +285,16 @@ SUFFIX(go_range_median_inter_nonconst) (DOUBLE *xs, int n, DOUBLE *res)
 }
 
 int
+SUFFIX(go_range_constant) (DOUBLE const *xs, int n)
+{
+	int i;
+	for (i = 1; i < n; i++)
+		if (xs[0] != xs[i])
+			return 0;
+	return 1;
+}
+
+int
 SUFFIX(go_range_increasing) (DOUBLE const *xs, int n)
 {
 	int i;
diff --git a/goffice/math/go-rangefunc.h b/goffice/math/go-rangefunc.h
index 006ccfb..4c9ff4c 100644
--- a/goffice/math/go-rangefunc.h
+++ b/goffice/math/go-rangefunc.h
@@ -17,6 +17,7 @@ int go_range_fractile_inter_sorted (double const *xs, int n, double *res, double
 int go_range_median_inter (double const *xs, int n, double *res);
 int go_range_median_inter_nonconst (double *xs, int n, double *res);
 int go_range_median_inter_sorted (double const *xs, int n, double *res);
+int go_range_constant (double const *xs, int n);
 int go_range_increasing (double const *xs, int n);
 int go_range_decreasing (double const *xs, int n);
 int go_range_vary_uniformly (double const *xs, int n);
@@ -36,6 +37,7 @@ int go_range_fractile_inter_sortedl (long double const *xs, int n, long double *
 int go_range_median_interl (long double const *xs, int n, long double *res);
 int go_range_median_inter_nonconstl (long double *xs, int n, long double *res);
 int go_range_median_inter_sortedl (long double const *xs, int n, long double *res);
+int go_range_constantl (long double const *xs, int n);
 int go_range_increasingl (long double const *xs, int n);
 int go_range_decreasingl (long double const *xs, int n);
 int go_range_vary_uniformlyl (long double const *xs, int n);
diff --git a/goffice/math/goffice-math.h b/goffice/math/goffice-math.h
index c307224..d0ae6f0 100644
--- a/goffice/math/goffice-math.h
+++ b/goffice/math/goffice-math.h
@@ -1,15 +1,16 @@
 #ifndef __GOFFICE_MATH_H
 #define __GOFFICE_MATH_H
 
+#include <goffice/math/go-accumulator.h>
 #include <goffice/math/go-complex.h>
 #include <goffice/math/go-cspline.h>
 #include <goffice/math/go-distribution.h>
 #include <goffice/math/go-fft.h>
 #include <goffice/math/go-math.h>
 #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>
+#include <goffice/math/go-rangefunc.h>
+#include <goffice/math/go-regression.h>
 
 #endif
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]