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 #include #include +#include namespace util { /////////////////////////////////////////////////////////////////////// @@ -741,7 +743,17 @@ namespace util { T sum (const K 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::infinity (), + ::util::distance (a, b), + "distance with an infinity is infinite" + ); + } return tap.status (); }