[goffice] range functions: improve precision when the input vector is (near-)constant.

[Morten Welinder <mortenw src gnome org>]

commit b639a1c505b868f6ab5a55d74d5690172ff01489
Author: Morten Welinder <terra gnome org>
Date:   Fri Sep 30 15:16:19 2011 -0400

    range functions: improve precision when the input vector is (near-)constant.

 ChangeLog                   |    3 ++
 NEWS                        |    1 +
 goffice/math/go-rangefunc.c |   53 +++++++++++++++++++++++++++++++++++++-----
 3 files changed, 50 insertions(+), 7 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 51f7630..26abe7f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,8 @@
 2011-09-30  Morten Welinder  <terra gnome org>
 
+	* goffice/math/go-rangefunc.c (sum_helper): Improve precision in
+	the all-constant case.  See #660564.
+
 	* goffice/utils/go-file.c (go_dirname_from_uri): Plug leak.
 
 2011-09-26  Jean Brefort  <jean brefort normalesup org>
diff --git a/NEWS b/NEWS
index 6c61d34..e89b065 100644
--- a/NEWS
+++ b/NEWS
@@ -31,6 +31,7 @@ Jean:
 Morten:
 	* Recognize scientific formats with longer exponents.
 	* Fix problem with sticky format colouring.
+	* Improve certain range functions in corner cases.  [#660564]
 
 --------------------------------------------------------------------------
 goffice 0.8.17:
diff --git a/goffice/math/go-rangefunc.c b/goffice/math/go-rangefunc.c
index 8558234..e62bb44 100644
--- a/goffice/math/go-rangefunc.c
+++ b/goffice/math/go-rangefunc.c
@@ -42,13 +42,41 @@
 
 /* ------------------------------------------------------------------------- */
 
+static int
+SUFFIX(identical_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)
+SUFFIX(sum_helper) (DOUBLE const *xs, int n, int *all_id)
 {
-	LDOUBLE sum = 0;
+	LDOUBLE sum;
 	int i;
 
-	for (i = 0; i < n; 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 ? 0 : i * (LDOUBLE)(xs[0]);
+
+	for (; i < n; i++)
 		sum += xs[i];
 
 	return sum;
@@ -60,7 +88,8 @@ SUFFIX(sum_helper) (DOUBLE const *xs, int n)
 int
 SUFFIX(go_range_sum) (DOUBLE const *xs, int n, DOUBLE *res)
 {
-	*res = SUFFIX(sum_helper) (xs, n);
+	int all_id;
+	*res = SUFFIX(sum_helper) (xs, n, &all_id);
 	return 0;
 }
 
@@ -89,10 +118,15 @@ 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) / n;
+	*res = SUFFIX(sum_helper) (xs, n, &all_id) / n;
+	if (all_id)
+		*res = xs[0];
+
 	return 0;
 }
 
@@ -154,8 +188,13 @@ SUFFIX(go_range_devsq) (DOUBLE const *xs, int n, DOUBLE *res)
 	LDOUBLE q = 0;
 
 	if (n > 0) {
-		int i;
-		LDOUBLE m = SUFFIX(sum_helper) (xs, n) / n;
+		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;