add a function to return the kahan sum of an array

For some reason the fit seems to have trouble with the kinetic energy.
Basically, it seems to "converge" even though when you run the minimization
again it finds a better minimum with a lower energy. I think this is likely due
to the fact that for muons the kinetic energy only really affects the range of
the muon and this is subject to error in the numerical integration.

I also thought that maybe it could be due to roundoff error in the likelihood
calculation, so I implemented the Kahan summation to try and reduce that. No
idea if it's actually improving things, but I should benchmark it later to see.

1 files changed, 20 insertions, 0 deletions

diff --git a/src/misc.c b/src/misc.c
index 2f02595..bff0ba9 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -2,6 +2,26 @@
 #include <math.h>
 #include <stdlib.h> /* for size_t */
 
+double kahan_sum(double *x, size_t n)
+{
+    /* Returns the sum of the elements of `x` using the Kahan summation algorithm.
+     *
+     * See https://en.wikipedia.org/wiki/Kahan_summation_algorithm. */
+    size_t i;
+    double sum, c, y, t;
+
+    sum = 0.0;
+    c = 0.0;
+    for (i = 0; i < n; i++) {
+        y = x[i] - c;
+        t = sum + y;
+        c = (t - sum) - y;
+        sum = t;
+    }
+
+    return sum;
+}
+
 double interp1d(double x, double *xp, double *yp, size_t n)
 {
     /* A fast interpolation routine which assumes that the values in `xp` are