diff --git a/vtkm/Math.h b/vtkm/Math.h index e9f762eaa..18e86932e 100644 --- a/vtkm/Math.h +++ b/vtkm/Math.h @@ -2701,12 +2701,6 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32 template inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d) { - // This is a bit awkward. clang-11 does not define FP_FAST_FMA, but still compiles this to the correct assembly. - // Windows, however, generates truly horrendous assembly from this, with no fmas, to the extent I assume it could - // contort itself into a performance bug. - // That said, MSVC converts even a*b - c*d into horrible assembly, so it may be a wash. - // You'd want to just use #ifdef FP_FAST_FMA, but then you'd lose the (correct) generated assembly on clang. - // See: https://stackoverflow.com/a/40765925/904050 T cd = c * d; T err = std::fma(-c, d, cd); T dop = std::fma(a, b, -cd); @@ -2719,7 +2713,7 @@ inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d) // If there are no real roots, both elements are NaNs. // The error should be at most 3 ulps. template -inline VTKM_EXEC_CONT vtkm::Pair QuadraticRoots(T a, T b, T c) +inline VTKM_EXEC_CONT vtkm::Vec QuadraticRoots(T a, T b, T c) { if (a == 0) { @@ -2728,19 +2722,19 @@ inline VTKM_EXEC_CONT vtkm::Pair QuadraticRoots(T a, T b, T c) if (c == 0) { // A degenerate case. All real numbers are roots; hopefully this arbitrary decision interacts gracefully with use. - return vtkm::Pair(0,0); + return vtkm::Vec(0,0); } else { - return vtkm::Pair(vtkm::Nan(), vtkm::Nan()); + return vtkm::Vec(vtkm::Nan(), vtkm::Nan()); } } - return vtkm::Pair(-c/b, -c/b); + return vtkm::Vec(-c/b, -c/b); } T delta = DifferenceOfProducts(b, b, 4*a, c); if (delta < 0) { - return vtkm::Pair(vtkm::Nan(), vtkm::Nan()); + return vtkm::Vec(vtkm::Nan(), vtkm::Nan()); } T q = -(b + vtkm::CopySign(vtkm::Sqrt(delta), b)) / 2; @@ -2748,9 +2742,9 @@ inline VTKM_EXEC_CONT vtkm::Pair QuadraticRoots(T a, T b, T c) T r1 = c / q; if (r0 < r1) { - return vtkm::Pair(r0, r1); + return vtkm::Vec(r0, r1); } - return vtkm::Pair(r1, r0); + return vtkm::Vec(r1, r0); } /// Bitwise operations diff --git a/vtkm/Math.h.in b/vtkm/Math.h.in index b8c5d39c6..48158890f 100644 --- a/vtkm/Math.h.in +++ b/vtkm/Math.h.in @@ -1303,12 +1303,6 @@ inline VTKM_EXEC_CONT vtkm::UInt64 FloatDistance(vtkm::Float32 x, vtkm::Float32 template inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d) { - // This is a bit awkward. clang-11 does not define FP_FAST_FMA, but still compiles this to the correct assembly. - // Windows, however, generates truly horrendous assembly from this, with no fmas, to the extent I assume it could - // contort itself into a performance bug. - // That said, MSVC converts even a*b - c*d into horrible assembly, so it may be a wash. - // You'd want to just use #ifdef FP_FAST_FMA, but then you'd lose the (correct) generated assembly on clang. - // See: https://stackoverflow.com/a/40765925/904050 T cd = c * d; T err = std::fma(-c, d, cd); T dop = std::fma(a, b, -cd); @@ -1321,7 +1315,7 @@ inline VTKM_EXEC_CONT T DifferenceOfProducts(T a, T b, T c, T d) // If there are no real roots, both elements are NaNs. // The error should be at most 3 ulps. template -inline VTKM_EXEC_CONT vtkm::Pair QuadraticRoots(T a, T b, T c) +inline VTKM_EXEC_CONT vtkm::Vec QuadraticRoots(T a, T b, T c) { if (a == 0) { @@ -1330,19 +1324,19 @@ inline VTKM_EXEC_CONT vtkm::Pair QuadraticRoots(T a, T b, T c) if (c == 0) { // A degenerate case. All real numbers are roots; hopefully this arbitrary decision interacts gracefully with use. - return vtkm::Pair(0,0); + return vtkm::Vec(0,0); } else { - return vtkm::Pair(vtkm::Nan(), vtkm::Nan()); + return vtkm::Vec(vtkm::Nan(), vtkm::Nan()); } } - return vtkm::Pair(-c/b, -c/b); + return vtkm::Vec(-c/b, -c/b); } T delta = DifferenceOfProducts(b, b, 4*a, c); if (delta < 0) { - return vtkm::Pair(vtkm::Nan(), vtkm::Nan()); + return vtkm::Vec(vtkm::Nan(), vtkm::Nan()); } T q = -(b + vtkm::CopySign(vtkm::Sqrt(delta), b)) / 2; @@ -1350,9 +1344,9 @@ inline VTKM_EXEC_CONT vtkm::Pair QuadraticRoots(T a, T b, T c) T r1 = c / q; if (r0 < r1) { - return vtkm::Pair(r0, r1); + return vtkm::Vec(r0, r1); } - return vtkm::Pair(r1, r0); + return vtkm::Vec(r1, r0); } /// Bitwise operations diff --git a/vtkm/testing/UnitTestMath.cxx b/vtkm/testing/UnitTestMath.cxx index 87f2813f7..35a01ddc1 100644 --- a/vtkm/testing/UnitTestMath.cxx +++ b/vtkm/testing/UnitTestMath.cxx @@ -916,7 +916,6 @@ struct ScalarVectorFieldTests : public vtkm::exec::FunctorBase vtkm::Float32 expected = 5.376600027084351f; vtkm::UInt64 dist = vtkm::FloatDistance(expected, computed); - std::cout << "Dist = " << dist << "\n"; VTKM_MATH_ASSERT( dist < 2, "Float distance for difference of products is which exceeds 1.5; this is in violation of a " @@ -934,17 +933,17 @@ void TestQuadraticRoots() const // (x-1)(x+1) = x² - 1: auto roots = vtkm::QuadraticRoots(1.0f, 0.0f, -1.0f); - vtkm::UInt64 dist = vtkm::FloatDistance(-1.0f, roots.first); + vtkm::UInt64 dist = vtkm::FloatDistance(-1.0f, roots[0]); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); - dist = vtkm::FloatDistance(1.0f, roots.second); + dist = vtkm::FloatDistance(1.0f, roots[1]); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); // No real roots: roots = vtkm::QuadraticRoots(1.0f, 0.0f, 1.0f); - VTKM_MATH_ASSERT(vtkm::IsNan(roots.first), + VTKM_MATH_ASSERT(vtkm::IsNan(roots[0]), "Roots should be Nan for a quadratic with complex roots."); - VTKM_MATH_ASSERT(vtkm::IsNan(roots.second), + VTKM_MATH_ASSERT(vtkm::IsNan(roots[1]), "Roots should be Nan for a quadratic with complex roots."); #ifdef FP_FAST_FMA @@ -952,18 +951,18 @@ void TestQuadraticRoots() const // x² + 200x - 0.000015 = 0 has roots // -200.000000075, 7.5e-8 roots = vtkm::QuadraticRoots(1.0f, 200.0f, -0.000015f); - dist = vtkm::FloatDistance(-200.000000075f, roots.first); + dist = vtkm::FloatDistance(-200.000000075f, roots[0]); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); - dist = vtkm::FloatDistance(7.5e-8f, roots.second); + dist = vtkm::FloatDistance(7.5e-8f, roots[1]); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); // Kahan's example: auto roots64 = vtkm::QuadraticRoots(94906265.625, 94906267.000, 94906268.375); - dist = vtkm::FloatDistance(1.0, roots64.first); + dist = vtkm::FloatDistance(1.0, roots64[0]); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); - dist = vtkm::FloatDistance(1.000000028975958, roots64.second); + dist = vtkm::FloatDistance(1.000000028975958, roots64[1]); VTKM_MATH_ASSERT(dist < 3, "Float distance for quadratic roots exceeds 3 ulps."); #endif }