From 046a369f55d27219fb78565ed9b5eaca72abf809 Mon Sep 17 00:00:00 2001
From: Danny Robson <danny@nerdcruft.net>
Date: Fri, 16 Jun 2017 17:38:55 +1000
Subject: [PATCH] coord: don't use kahan summation for sum(coord)

Our kahan summation algorithm has issues with infinities so we just
avoid using it in coordinates where speed and simplicity are more
important.
---
 coord/ops.hpp  | 14 +++++++++++++-
 maths.hpp      |  6 ++++++
 test/coord.cpp | 10 ++++++++++
 3 files changed, 29 insertions(+), 1 deletion(-)

diff --git a/coord/ops.hpp b/coord/ops.hpp
index e66a945c..765cb77f 100644
--- a/coord/ops.hpp
+++ b/coord/ops.hpp
@@ -24,8 +24,10 @@
 #include "../preprocessor.hpp"
 #include "../types/bits.hpp"
 
+#include <algorithm>
 #include <cmath>
 #include <cstdlib>
+#include <iterator>
 
 namespace util {
     ///////////////////////////////////////////////////////////////////////
@@ -741,7 +743,17 @@ namespace util {
     T
     sum (const K<S,T> k)
     {
-        return sum (std::cbegin (k), std::cend (k));
+        // DO NOT USE util::sum(begin, end) from maths.hpp
+        //
+        // It would be nice to use kahan summation from maths.hpp but speed
+        // and simplicity is more important for these fixed sized
+        // coordinates. Infinities tend to crop up using these classes and
+        // they cause a few headaches in the kahan code.
+        //
+        // So, if the user wants kahan summation they can request it
+        // explicitly.
+
+        return std::accumulate (std::cbegin (k), std::cend (k), T{0});
     }
 
 
diff --git a/maths.hpp b/maths.hpp
index 747f9f7f..9e8e9f35 100644
--- a/maths.hpp
+++ b/maths.hpp
@@ -549,8 +549,14 @@ namespace util {
         T c = 0;
 
         for (auto cursor = first; cursor != last; ++cursor) {
+            // Infinities are handled poorly in this implementation. We tend
+            // to produce NaNs because of the subtraction where we compute
+            // `c'. For the time being just panic in this scenario.
+            assert(!std::isinf (*cursor));
+
             T y = *cursor - c;
             T t = sum + y;
+
             c = (t - sum) - y;
             sum = t;
         }
diff --git a/test/coord.cpp b/test/coord.cpp
index d9c4d9b0..dfdfdb12 100644
--- a/test/coord.cpp
+++ b/test/coord.cpp
@@ -42,6 +42,16 @@ main (void)
         tap.expect (x == p.x && y == p.y, "structured bindings extract correct data");
     }
 
+    {
+        util::point3f a { 103, 0, 14 };
+        util::point3f b { 104, INFINITY, 15 };
+
+        tap.expect_eq (
+            std::numeric_limits<float>::infinity (),
+            ::util::distance (a, b),
+            "distance with an infinity is infinite"
+        );
+    }
 
     return tap.status ();
 }