[gnumeric] MMULT: Use GnmAccumulator to avoid certain rounding errors.



commit 57ab2325bef16349b419324224618608da83b108
Author: Morten Welinder <terra gnome org>
Date:   Wed Jul 25 18:26:04 2012 -0400

    MMULT: Use GnmAccumulator to avoid certain rounding errors.

 ChangeLog      |    6 ++++--
 src/mathfunc.c |   18 +++++++++++++-----
 src/numbers.h  |    2 ++
 3 files changed, 19 insertions(+), 7 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 41238cb..99428de 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,10 +1,12 @@
+2012-07-25  Morten Welinder  <terra gnome org>
+
+	* src/mathfunc.c (mmult): Use GnmAccumulator for extra accuracy.
+
 2012-07-25  Andreas J. Guelzow <aguelzow pyrshep ca>
 
 	* src/wbc-gtk-edit.c (wbcg_edit_finish): reset the expr_txt pointer
 	when changing the txt pointer. Fixes #680548
 
-2012-07-24  Morten Welinder  <terra gnome org>
-
 	* src/dependent.c (link_unlink_expr_dep): Fix problem from
 	earlier cleanup.
 
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 1e80b62..c6beb7a 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -6447,17 +6447,25 @@ void
 mmult (gnm_float *A, gnm_float *B, int cols_a, int rows_a, int cols_b,
        gnm_float *product)
 {
-	gnm_float tmp;
 	int	c, r, i;
+	void *state = gnm_accumulator_start ();
+	GnmAccumulator *acc = gnm_accumulator_new ();
 
 	for (c = 0; c < cols_b; ++c) {
 		for (r = 0; r < rows_a; ++r) {
-			tmp = 0;
-			for (i = 0; i < cols_a; ++i)
-				tmp += A[r + i * rows_a] * B[i + c * cols_a];
-			product[r + c * rows_a] = tmp;
+			go_accumulator_clear (acc);
+			for (i = 0; i < cols_a; ++i) {
+				GnmQuad p;
+				gnm_quad_mul12 (&p,
+						A[r + i * rows_a],
+						B[i + c * cols_a]);
+				gnm_accumulator_add_quad (acc, &p);
+			}
+			product[r + c * rows_a] = gnm_accumulator_value (acc);
 		}
 	}
+	gnm_accumulator_free (acc);
+	gnm_accumulator_end (state);
 }
 
 /***************************************************************************/
diff --git a/src/numbers.h b/src/numbers.h
index 112f339..58c76e9 100644
--- a/src/numbers.h
+++ b/src/numbers.h
@@ -123,6 +123,7 @@ gnm_float gnm_yn (int n, gnm_float x);
 #define GnmQuad GOQuadl
 #define gnm_quad_init go_quad_initl
 #define gnm_quad_mul go_quad_mull
+#define gnm_quad_mul12 go_quad_mul12l
 #define GnmAccumulator GOAccumulatorl
 #define gnm_accumulator_start go_accumulator_startl
 #define gnm_accumulator_end go_accumulator_endl
@@ -205,6 +206,7 @@ typedef double gnm_float;
 #define GnmQuad GOQuad
 #define gnm_quad_init go_quad_init
 #define gnm_quad_mul go_quad_mul
+#define gnm_quad_mul12 go_quad_mul12
 #define GnmAccumulator GOAccumulator
 #define gnm_accumulator_start go_accumulator_start
 #define gnm_accumulator_end go_accumulator_end



[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]