diff --git a/singularity-eos/base/fast-math/logs.hpp b/singularity-eos/base/fast-math/logs.hpp index 09dc2be059..7eb02bef62 100644 --- a/singularity-eos/base/fast-math/logs.hpp +++ b/singularity-eos/base/fast-math/logs.hpp @@ -67,11 +67,29 @@ PORTABLE_FORCEINLINE_FUNCTION auto sgn(const T &x) { return (T(0) < x) - (x < T(0)); } +// TODO(JMM): Add templating and concepts to enforce bit width +// equivalence in input/output types. PORTABLE_FORCEINLINE_FUNCTION auto as_int(double f) { return *reinterpret_cast(&f); } PORTABLE_FORCEINLINE_FUNCTION auto as_double(std::int64_t i) { return *reinterpret_cast(&i); } +// Magic numbers constexpr because C++ doesn't constexpr reinterpret casts +// These numbers will be different for different precisions and endianness. +namespace FP64LE { // 64 bit, little endian +constexpr std::int64_t one = 1; +// as_int(1.0) == 2^62 - 2^52 +constexpr std::int64_t one_as_int = (one << 62) - (one << 52); +// 1./static_cast(as_int(2.0) - as_int(1.0)) == 2^-52 +constexpr double scale_down = 2.22044604925031e-16; +// as_int(2.0) - as_int(1.0) = 2^52, but note the type +constexpr double scale_up = (one << 52); +// 2^52 - 1 +constexpr std::int64_t mantissa_mask = (one << 52) - one; +// 2^26 - 1 +constexpr std::int64_t low_mask = (one << 26) - 1; +} // namespace FP64LE + // First order interpolation based NQTs // ---------------------------------------------------------------------- // Reference implementations, however the integer cast implementation @@ -96,24 +114,14 @@ double pow2_o1_portable(const double x) { // Integer aliased versions PORTABLE_FORCEINLINE_FUNCTION double lg_o1_aliased(const double x) { + using namespace FP64LE; PORTABLE_REQUIRE(x > 0, "log divergent for x <= 0"); - // Magic numbers constexpr because C++ doesn't constexpr reinterpret casts - // these are floating point numbers as reinterpreted as integers. - // as_int(1.0) - constexpr std::int64_t one_as_int = 4607182418800017408; - // 1./static_cast(as_int(2.0) - as_int(1.0)) - constexpr double scale_down = 2.22044604925031e-16; return static_cast(as_int(x) - one_as_int) * scale_down; } PORTABLE_FORCEINLINE_FUNCTION double pow2_o1_aliased(const double x) { - // Magic numbers constexpr because C++ doesn't constexpr reinterpret casts - // these are floating point numbers as reinterpreted as integers. - // as_int(1.0) - constexpr std::int64_t one_as_int = 4607182418800017408; - // as_int(2.0) - as_int(1.0) - constexpr double scale_up = 4503599627370496; + using namespace FP64LE; return as_double(static_cast(x * scale_up) + one_as_int); } // ---------------------------------------------------------------------- @@ -145,16 +153,8 @@ double pow2_o2_portable(const double x) { // Integer aliased/bithacked versions PORTABLE_FORCEINLINE_FUNCTION double lg_o2_aliased(const double x) { + using namespace FP64LE; PORTABLE_REQUIRE(x > 0, "log divergent for x <= 0"); - // as_int(1.0) == 2^62 - 2^52 - constexpr std::int64_t one_as_int = 4607182418800017408; - // 1/(as_int(2.0) - as_int(1.0)) == 2^-52 - constexpr double scale_down = 2.220446049250313e-16; - // 2^52 - 1 - constexpr std::int64_t mantissa_mask = 4503599627370495; - // 2^26 - 1 - constexpr std::int64_t low_mask = 67108863; - const std::int64_t x_as_int = as_int(x) - one_as_int; const std::int64_t frac_as_int = x_as_int & mantissa_mask; const std::int64_t frac_high = frac_as_int >> 26; @@ -167,15 +167,10 @@ double lg_o2_aliased(const double x) { PORTABLE_FORCEINLINE_FUNCTION double pow2_o2_aliased(const double x) { - // as_int(1.0) == 2^62 - 2^52 - constexpr std::int64_t one_as_int = 4607182418800017408; - // as_int(2.0) - as_int(1.0) == 2^52 - constexpr double scale_up = 4503599627370496; - constexpr std::int64_t mantissa_mask = 4503599627370495; // 2^52 - 1 - constexpr std::int64_t a = 9007199254740992; // 2 * 2^52 - constexpr double b = 67108864; // 2^26 - constexpr std::int64_t c = 18014398509481984; // 4 * 2^52 - + using namespace FP64LE; + constexpr std::int64_t a = 9007199254740992; // 2 * 2^52 + constexpr double b = 67108864; // 2^26 + constexpr std::int64_t c = 18014398509481984; // 4 * 2^52 const std::int64_t x_as_int = static_cast(x * scale_up); const std::int64_t frac_as_int = x_as_int & mantissa_mask; const std::int64_t frac_sqrt = diff --git a/test/test_math_utils.cpp b/test/test_math_utils.cpp index feff8316bc..1dee262ef6 100644 --- a/test/test_math_utils.cpp +++ b/test/test_math_utils.cpp @@ -41,6 +41,18 @@ SCENARIO("Test that we can either throw an error on host or do nothing on device REQUIRE_MAYBE_THROWS(PORTABLE_ALWAYS_THROW_OR_ABORT("Error message")); } +#ifndef SINGULARITY_USE_TRUE_LOG_GRIDDING +#ifndef SINGULARITY_NQT_PORTABLE +SCENARIO("Test that the fast math magic numbers are all correct", "[FastMath]") { + REQUIRE(singularity::FastMath::FP64LE::one_as_int == 4607182418800017408); + REQUIRE(singularity::FastMath::FP64LE::scale_down == 2.22044604925031e-16); + REQUIRE(singularity::FastMath::FP64LE::scale_up == 4503599627370496.0); + REQUIRE(singularity::FastMath::FP64LE::mantissa_mask == 4503599627370495); + REQUIRE(singularity::FastMath::FP64LE::low_mask == 67108863); +} +#endif +#endif + SCENARIO("Test that fast logs are invertible and run on device", "[FastMath]") { GIVEN("A set of values to invert over a large dynamic range") { constexpr Real LXMIN = -20; @@ -64,7 +76,7 @@ SCENARIO("Test that fast logs are invertible and run on device", "[FastMath]") { portableFor( "try out the fast math", 0, NX, PORTABLE_LAMBDA(const int i) { constexpr Real machine_eps = std::numeric_limits::epsilon(); - constexpr Real acceptable_err = 100 * machine_eps; + constexpr Real acceptable_err = 1000 * machine_eps; const Real lx = singularity::FastMath::log10(x[i]); const Real elx = singularity::FastMath::pow10(lx); const Real rel_err = 2.0 * std::abs(x[i] - elx) /