From 470058c328c5772b01003fb3accc6892f02077d0 Mon Sep 17 00:00:00 2001 From: Nick Date: Wed, 24 Jun 2020 10:52:59 -0400 Subject: [PATCH] Unfoobar the float_distance MR. --- vtkm/Math.h | 101 +++++++++++++++ vtkm/Math.h.in | 101 +++++++++++++++ vtkm/testing/TestingMath.h | 259 +++++++++++++++++++++++++++++++++++++ 3 files changed, 461 insertions(+) diff --git a/vtkm/Math.h b/vtkm/Math.h index 8c7ddc8f1..7525d2dc9 100644 --- a/vtkm/Math.h +++ b/vtkm/Math.h @@ -17,8 +17,10 @@ #include #include +#include // must be found with or without CUDA. #ifndef VTKM_CUDA #include +#include #include #include #include @@ -2584,6 +2586,105 @@ inline VTKM_EXEC_CONT vtkm::Float64 Ldexp(vtkm::Float64 x, vtkm::Int32 exponent) #endif } +// See: https://randomascii.wordpress.com/2012/01/23/stupid-float-tricks-2/ for why this works. +inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float64 x, vtkm::Float64 y) +{ + static_assert(sizeof(vtkm::Float64) == sizeof(vtkm::UInt64), "vtkm::Float64 is incorrect size."); + static_assert(std::numeric_limits::has_denorm == std::denorm_present, "FloatDistance presumes the floating-point type has subnormal numbers."); + + if (!vtkm::IsFinite(x) || !vtkm::IsFinite(y)) { + return 0xFFFFFFFFFFFFFFFFL; + } + + // Signed zero is the sworn enemy of this process. + if (y == 0) { + y = vtkm::Abs(y); + } + if (x == 0) { + x = vtkm::Abs(x); + } + + if ( (x < 0 && y >= 0) || (x >= 0 && y < 0) ) + { + vtkm::UInt64 dx, dy; + if (x < 0) { + dy = FloatDistance(0.0, y); + dx = FloatDistance(0.0, -x); + } + else { + dy = FloatDistance(0.0, -y); + dx = FloatDistance(0.0, x); + } + + return dx + dy; + } + + if (x < 0 && y < 0) { + return FloatDistance(-x, -y); + } + + // Note that: + // int64_t xi = *reinterpret_cast(&x); + // int64_t yi = *reinterpret_cast(&y); + // also works, but generates warnings. + // Good option to have if we get compile errors off memcpy or don't want to #include though. + // At least on gcc, both versions generate the same assembly. + vtkm::UInt64 xi; + vtkm::UInt64 yi; + memcpy(&xi, &x, sizeof(vtkm::UInt64)); + memcpy(&yi, &y, sizeof(vtkm::UInt64)); + if (yi > xi) { + return yi - xi; + } + return xi - yi; +} + +inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32 y) +{ + static_assert(sizeof(vtkm::Float32) == sizeof(vtkm::Int32), "vtkm::Float32 is incorrect size."); + static_assert(std::numeric_limits::has_denorm == std::denorm_present, "FloatDistance presumes the floating-point type has subnormal numbers."); + + if (!vtkm::IsFinite(x) || !vtkm::IsFinite(y)) { + return 0xFFFFFFFFFFFFFFFFL; + } + + if (y == 0) { + y = vtkm::Abs(y); + } + if (x == 0) { + x = vtkm::Abs(x); + } + + if ( (x < 0 && y >= 0) || (x >= 0 && y < 0) ) + { + vtkm::UInt64 dx, dy; + if (x < 0) { + dy = FloatDistance(0.0f, y); + dx = FloatDistance(0.0f, -x); + } + else { + dy = FloatDistance(0.0f, -y); + dx = FloatDistance(0.0f, x); + } + return dx + dy; + } + + if (x < 0 && y < 0) { + return FloatDistance(-x, -y); + } + + vtkm::UInt32 xi_32; + vtkm::UInt32 yi_32; + memcpy(&xi_32, &x, sizeof(vtkm::UInt32)); + memcpy(&yi_32, &y, sizeof(vtkm::UInt32)); + vtkm::UInt64 xi = xi_32; + vtkm::UInt64 yi = yi_32; + if (yi > xi) { + return yi - xi; + } + return xi - yi; +} + /// Bitwise operations /// diff --git a/vtkm/Math.h.in b/vtkm/Math.h.in index 12a64ce7e..3fb1f30de 100644 --- a/vtkm/Math.h.in +++ b/vtkm/Math.h.in @@ -29,8 +29,10 @@ $# Ignore the following comment. It is meant for the generated file. #include #include +#include // must be found with or without CUDA. #ifndef VTKM_CUDA #include +#include #include #include #include @@ -1186,6 +1188,105 @@ inline VTKM_EXEC_CONT vtkm::Float64 Ldexp(vtkm::Float64 x, vtkm::Int32 exponent) #endif } +// See: https://randomascii.wordpress.com/2012/01/23/stupid-float-tricks-2/ for why this works. +inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float64 x, vtkm::Float64 y) +{ + static_assert(sizeof(vtkm::Float64) == sizeof(vtkm::UInt64), "vtkm::Float64 is incorrect size."); + static_assert(std::numeric_limits::has_denorm == std::denorm_present, "FloatDistance presumes the floating-point type has subnormal numbers."); + + if (!vtkm::IsFinite(x) || !vtkm::IsFinite(y)) { + return 0xFFFFFFFFFFFFFFFFL; + } + + // Signed zero is the sworn enemy of this process. + if (y == 0) { + y = vtkm::Abs(y); + } + if (x == 0) { + x = vtkm::Abs(x); + } + + if ( (x < 0 && y >= 0) || (x >= 0 && y < 0) ) + { + vtkm::UInt64 dx, dy; + if (x < 0) { + dy = FloatDistance(0.0, y); + dx = FloatDistance(0.0, -x); + } + else { + dy = FloatDistance(0.0, -y); + dx = FloatDistance(0.0, x); + } + + return dx + dy; + } + + if (x < 0 && y < 0) { + return FloatDistance(-x, -y); + } + + // Note that: + // int64_t xi = *reinterpret_cast(&x); + // int64_t yi = *reinterpret_cast(&y); + // also works, but generates warnings. + // Good option to have if we get compile errors off memcpy or don't want to #include though. + // At least on gcc, both versions generate the same assembly. + vtkm::UInt64 xi; + vtkm::UInt64 yi; + memcpy(&xi, &x, sizeof(vtkm::UInt64)); + memcpy(&yi, &y, sizeof(vtkm::UInt64)); + if (yi > xi) { + return yi - xi; + } + return xi - yi; +} + +inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32 y) +{ + static_assert(sizeof(vtkm::Float32) == sizeof(vtkm::Int32), "vtkm::Float32 is incorrect size."); + static_assert(std::numeric_limits::has_denorm == std::denorm_present, "FloatDistance presumes the floating-point type has subnormal numbers."); + + if (!vtkm::IsFinite(x) || !vtkm::IsFinite(y)) { + return 0xFFFFFFFFFFFFFFFFL; + } + + if (y == 0) { + y = vtkm::Abs(y); + } + if (x == 0) { + x = vtkm::Abs(x); + } + + if ( (x < 0 && y >= 0) || (x >= 0 && y < 0) ) + { + vtkm::UInt64 dx, dy; + if (x < 0) { + dy = FloatDistance(0.0f, y); + dx = FloatDistance(0.0f, -x); + } + else { + dy = FloatDistance(0.0f, -y); + dx = FloatDistance(0.0f, x); + } + return dx + dy; + } + + if (x < 0 && y < 0) { + return FloatDistance(-x, -y); + } + + vtkm::UInt32 xi_32; + vtkm::UInt32 yi_32; + memcpy(&xi_32, &x, sizeof(vtkm::UInt32)); + memcpy(&yi_32, &y, sizeof(vtkm::UInt32)); + vtkm::UInt64 xi = xi_32; + vtkm::UInt64 yi = yi_32; + if (yi > xi) { + return yi - xi; + } + return xi - yi; +} + /// Bitwise operations /// diff --git a/vtkm/testing/TestingMath.h b/vtkm/testing/TestingMath.h index 864980d57..8c04e1384 100644 --- a/vtkm/testing/TestingMath.h +++ b/vtkm/testing/TestingMath.h @@ -646,6 +646,264 @@ struct ScalarVectorFieldTests : public vtkm::exec::FunctorBase "CopySign failed."); } + VTKM_EXEC + void TestFloatDistance() const + { + { + vtkm::UInt64 dist = vtkm::FloatDistance(1.0, 1.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from 1.0 to 1.0 is not zero."); + + dist = vtkm::FloatDistance(-1.0, -1.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from -1.0 to -1.0 is not zero."); + + dist = vtkm::FloatDistance(0.0, 0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from 0.0 to 0.0 is not zero."); + + // Check nan: + dist = vtkm::FloatDistance(std::numeric_limits::quiet_NaN(), 1.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0xFFFFFFFFFFFFFFFFL), dist), + "Float distance to a Nan is not the documented value."); + + dist = vtkm::FloatDistance(1.0, std::numeric_limits::quiet_NaN()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0xFFFFFFFFFFFFFFFFL), dist), + "Float distance to a Nan is not the documented value."); + + // Check infinity: + dist = vtkm::FloatDistance(std::numeric_limits::infinity(), 1.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0xFFFFFFFFFFFFFFFFL), dist), + "Float distance to infinity is not the documented value."); + + dist = vtkm::FloatDistance(1.0, std::numeric_limits::infinity()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0xFFFFFFFFFFFFFFFFL), dist), + "Float distance to infinity is not the documented value."); + + // Check saturation: + dist = vtkm::FloatDistance(std::numeric_limits::lowest(), + std::numeric_limits::max()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(18437736874454810622uL), dist), + "Float distance from lowest to max is incorrect."); + + dist = vtkm::FloatDistance(std::numeric_limits::max(), + std::numeric_limits::lowest()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(18437736874454810622uL), dist), + "Float distance from max to lowest is incorrect."); + + // Check symmetry: + dist = vtkm::FloatDistance(-2.0, -1.0); + vtkm::UInt64 dist2 = vtkm::FloatDistance(-1.0, -2.0); + VTKM_MATH_ASSERT(test_equal(dist2, dist), "Symmetry of negative numbers does not hold."); + + dist = vtkm::FloatDistance(1.0, 2.0); + dist2 = vtkm::FloatDistance(2.0, 1.0); + VTKM_MATH_ASSERT(test_equal(dist2, dist), "Float distance 1->2 != float distance 2->1."); + + // Check symmetry of bound which includes zero: + dist = vtkm::FloatDistance(-0.25, 0.25); + dist2 = vtkm::FloatDistance(0.25, -0.25); + VTKM_MATH_ASSERT(test_equal(dist2, dist), + "Symmetry is violated over a bound which contains zero."); + + // Check correctness: + dist = vtkm::FloatDistance(1.0, 1.0 + std::numeric_limits::epsilon()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), + "Float distance from 1 to 1 + eps is not = 1."); + dist = vtkm::FloatDistance(1.0 + std::numeric_limits::epsilon(), 1.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated"); + + dist = vtkm::FloatDistance(1.0, 1.0 + 2 * std::numeric_limits::epsilon()); + VTKM_MATH_ASSERT(test_equal(vtkm::Int64(2), dist), + "Float distance from 1 to 1 + 2eps is not 2."); + dist = vtkm::FloatDistance(1.0 + 2 * std::numeric_limits::epsilon(), 1.0); + VTKM_MATH_ASSERT(test_equal(vtkm::Int64(2), dist), "Symmetry is violated."); + + // Now test x = y: + vtkm::Float64 x = -1; + for (int i = 0; i < 500; ++i) + { + dist = vtkm::FloatDistance(x, x); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from x to x is not zero."); + x += 0.01; + } + // Test zero: + dist = vtkm::FloatDistance(0.0, 0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from zero to zero is not zero."); + // Test signed zero: + dist = vtkm::FloatDistance(0.0, -0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from 0.0 to -0.0 is not zero."); + + dist = vtkm::FloatDistance(-0.0, 0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from -0.0 to 0.0 is not zero."); + + dist = vtkm::FloatDistance(-0.0, -0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from -0.0 to 0.0 is not zero."); + + // Negative to negative zero: + dist = vtkm::FloatDistance(-std::numeric_limits::denorm_min(), -0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Negative to zero incorrect."); + // And symmetry: + dist = vtkm::FloatDistance(-0.0, -std::numeric_limits::denorm_min()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated."); + + // Negative to positive zero: + dist = vtkm::FloatDistance(-std::numeric_limits::denorm_min(), 0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), + "Negative to positive zero is incorrect."); + // And symmetry: + dist = vtkm::FloatDistance(0.0, -std::numeric_limits::denorm_min()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated."); + + // Positive to zero: + dist = vtkm::FloatDistance(std::numeric_limits::denorm_min(), 0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Positive to zero is incorrect."); + // And symmetry: + dist = vtkm::FloatDistance(0.0, std::numeric_limits::denorm_min()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated"); + + // Positive to negative zero: + dist = vtkm::FloatDistance(std::numeric_limits::denorm_min(), -0.0); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), + "Positive to negative zero is incorrect."); + // And symmetry: + dist = vtkm::FloatDistance(-0.0, std::numeric_limits::denorm_min()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated."); + } + + // I would try to just template these, but in fact the double precision version has to saturate, + // whereas the float version has sufficient range. + { + vtkm::UInt64 dist = vtkm::FloatDistance(1.0f, 1.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from 1.0 to 1.0 is not zero."); + + dist = vtkm::FloatDistance(-1.0f, -1.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from -1.0 to -1.0 is not zero."); + + dist = vtkm::FloatDistance(0.0f, 0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from 0.0 to 0.0 is not zero."); + + // Check nan: + dist = vtkm::FloatDistance(std::numeric_limits::quiet_NaN(), 1.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0xFFFFFFFFFFFFFFFFL), dist), + "Float distance to a Nan is not the documented value."); + + dist = vtkm::FloatDistance(1.0f, std::numeric_limits::quiet_NaN()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0xFFFFFFFFFFFFFFFFL), dist), + "Float distance to a Nan is not the documented value."); + + // Check infinity: + dist = vtkm::FloatDistance(std::numeric_limits::infinity(), 1.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0xFFFFFFFFFFFFFFFFL), dist), + "Float distance to infinity is not the documented value."); + + dist = vtkm::FloatDistance(1.0f, std::numeric_limits::infinity()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0xFFFFFFFFFFFFFFFFL), dist), + "Float distance to infinity is not the documented value."); + + // Check saturation: + dist = vtkm::FloatDistance(std::numeric_limits::lowest(), + std::numeric_limits::max()); + VTKM_MATH_ASSERT(dist > 0, "Float distance is negative."); + + dist = vtkm::FloatDistance(std::numeric_limits::max(), + std::numeric_limits::lowest()); + VTKM_MATH_ASSERT(dist > 0, "Float distance is negative."); + + // Check symmetry: + dist = vtkm::FloatDistance(-2.0f, -1.0f); + vtkm::UInt64 dist2 = vtkm::FloatDistance(-1.0f, -2.0f); + VTKM_MATH_ASSERT(test_equal(dist2, dist), "Symmetry of negative numbers does not hold."); + + dist = vtkm::FloatDistance(1.0f, 2.0f); + dist2 = vtkm::FloatDistance(2.0f, 1.0f); + VTKM_MATH_ASSERT(test_equal(dist2, dist), "Float distance 1->2 != float distance 2->1."); + + // Check symmetry of bound which includes zero: + dist = vtkm::FloatDistance(-0.25f, 0.25f); + dist2 = vtkm::FloatDistance(0.25f, -0.25f); + VTKM_MATH_ASSERT(test_equal(dist2, dist), + "Symmetry is violated over a bound which contains zero."); + + // Check correctness: + dist = vtkm::FloatDistance(1.0f, 1.0f + std::numeric_limits::epsilon()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), + "Float distance from 1 to 1 + eps is not = 1."); + dist = vtkm::FloatDistance(1.0f + std::numeric_limits::epsilon(), 1.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated"); + + dist = vtkm::FloatDistance(1.0f, 1.0f + 2 * std::numeric_limits::epsilon()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(2), dist), + "Float distance from 1 to 1 + 2eps is not 2."); + dist = vtkm::FloatDistance(1.0f + 2 * std::numeric_limits::epsilon(), 1.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(2), dist), "Symmetry is violated."); + + // Now test x = y: + vtkm::Float32 x = -1; + for (int i = 0; i < 500; ++i) + { + dist = vtkm::FloatDistance(x, x); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from x to x is not zero."); + x += 0.01f; + } + // Test zero: + dist = vtkm::FloatDistance(0.0f, 0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from zero to zero is not zero."); + // Test signed zero: + dist = vtkm::FloatDistance(0.0f, -0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from 0.0 to -0.0 is not zero."); + + dist = vtkm::FloatDistance(-0.0f, 0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from -0.0 to 0.0 is not zero."); + + dist = vtkm::FloatDistance(-0.0f, -0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(0), dist), + "Float distance from -0.0 to 0.0 is not zero."); + + // Negative to negative zero: + dist = vtkm::FloatDistance(-std::numeric_limits::denorm_min(), -0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Negative to zero incorrect."); + // And symmetry: + dist = vtkm::FloatDistance(-0.0f, -std::numeric_limits::denorm_min()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated."); + + // Negative to positive zero: + dist = vtkm::FloatDistance(-std::numeric_limits::denorm_min(), 0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), + "Negative to positive zero is incorrect."); + // And symmetry: + dist = vtkm::FloatDistance(0.0f, -std::numeric_limits::denorm_min()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated."); + + // Positive to zero: + dist = vtkm::FloatDistance(std::numeric_limits::denorm_min(), 0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Positive to zero is incorrect."); + // And symmetry: + dist = vtkm::FloatDistance(0.0f, std::numeric_limits::denorm_min()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated"); + + // Positive to negative zero: + dist = vtkm::FloatDistance(std::numeric_limits::denorm_min(), -0.0f); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), + "Positive to negative zero is incorrect."); + // And symmetry: + dist = vtkm::FloatDistance(-0.0f, std::numeric_limits::denorm_min()); + VTKM_MATH_ASSERT(test_equal(vtkm::UInt64(1), dist), "Symmetry is violated."); + } + } + VTKM_EXEC void operator()(vtkm::Id) const { @@ -663,6 +921,7 @@ struct ScalarVectorFieldTests : public vtkm::exec::FunctorBase this->TestLog10(); this->TestLog1P(); this->TestCopySign(); + this->TestFloatDistance(); } };