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, 29 insertions, 1 deletions

diff --git a/src/test.c b/src/test.c
index 2ea73d1..fb788ba 100644
--- a/src/test.c
+++ b/src/test.c
@@ -440,6 +440,33 @@ err:
     return 1;
 }
 
+int test_kahan_sum(char *err)
+{
+    /* Tests the kahan_sum() function. */
+    size_t i;
+    double x[100], sum, expected;
+
+    init_genrand(0);
+
+    expected = 0.0;
+    for (i = 0; i < sizeof(x)/sizeof(x[0]); i++) {
+        x[i] = genrand_real2();
+        expected += x[i];
+    }
+
+    sum = kahan_sum(x,sizeof(x)/sizeof(x[0]));
+
+    if (!isclose(sum, expected, 1e-9, 1e-9)) {
+        sprintf(err, "kahan_sum returned %.5g, but expected %.5g", sum, expected);
+        goto err;
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -453,7 +480,8 @@ struct tests {
     {test_logsumexp, "test_logsumexp"},
     {test_sno_charge_integral, "test_sno_charge_integral"},
     {test_path, "test_path"},
-    {test_interp1d, "test_interp1d"}
+    {test_interp1d, "test_interp1d"},
+    {test_kahan_sum, "test_kahan_sum"},
 };
 
 int main(int argc, char **argv)