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.
This commit is contained in:
Danny Robson 2017-06-16 17:38:55 +10:00
parent 061088467f
commit 046a369f55
3 changed files with 29 additions and 1 deletions

View File

@ -24,8 +24,10 @@
#include "../preprocessor.hpp" #include "../preprocessor.hpp"
#include "../types/bits.hpp" #include "../types/bits.hpp"
#include <algorithm>
#include <cmath> #include <cmath>
#include <cstdlib> #include <cstdlib>
#include <iterator>
namespace util { namespace util {
/////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////
@ -741,7 +743,17 @@ namespace util {
T T
sum (const K<S,T> k) 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});
} }

View File

@ -549,8 +549,14 @@ namespace util {
T c = 0; T c = 0;
for (auto cursor = first; cursor != last; ++cursor) { 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 y = *cursor - c;
T t = sum + y; T t = sum + y;
c = (t - sum) - y; c = (t - sum) - y;
sum = t; sum = t;
} }

View File

@ -42,6 +42,16 @@ main (void)
tap.expect (x == p.x && y == p.y, "structured bindings extract correct data"); 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 (); return tap.status ();
} }