diff --git a/source/tools/math.h b/source/tools/math.h index 8fe33467e2..69dff51d85 100644 --- a/source/tools/math.h +++ b/source/tools/math.h @@ -214,13 +214,15 @@ namespace emp { } } + namespace internal { /// A fast (O(log p)) integral-power command. template - static constexpr type_if Pow(T base, T p) { + static constexpr T PowIntImpl(T base, T p) { if (p <= 0) return 1; - if (p & 1) return base * Pow(base, p-1); // Odd exponent: strip one mulitple off and recurse. - return Square( Pow(base,p/2) ); // Even exponent: calc for half and square result. + if (p & 1) return base * PowIntImpl(base, p-1); // Odd exponent: strip one mulitple off and recurse. + return emp::Square( PowIntImpl(base,p/2) ); // Even exponent: calc for half and square result. } + } // namespace internal /// A fast 2^x command. static constexpr double Pow2(double exp) { @@ -234,11 +236,41 @@ namespace emp { return exp < 1 ? 1 : (base * IntPow(base, exp-1)); } + namespace internal { /// A fast method for calculating exponents on doubles. - static constexpr double Pow(double base, double exp) { + static constexpr double PowDoubleImpl(double base, double exp) { // Normally, convert to a base of 2 and then use Pow2. // If base is negative, we don't want to deal with imaginary numbers, so use IntPow. - return (base > 0) ? Pow2(Log2(base) * exp) : IntPow(base,exp); + return (base > 0) + ? emp::Pow2(emp::Log2(base) * exp) + : emp::IntPow(base,exp); + } + + // adapted from https://stackoverflow.com/a/30836042 + // prevents argument from being used for type deduction + template struct identity { typedef T type; }; + + } // namespace internal + + + /// A fast method for calculating exponents on doubles or integral types. + /// Uses if constexpr to work around compiler bug in Emscripten (issue #296). + template + static constexpr T Pow( + T base, typename internal::identity::type exp + ) { + // TODO cpp20, C++20 replace with std::is_constant_evaluated + // adapted from https://stackoverflow.com/a/62610143 + // exclude clang versions with compiler bug https://reviews.llvm.org/D35190 + #if defined(__clang__) && __clang_major__>=9 || defined(__GNUC__) && !defined(__clang__) + // if base is not known at compile time, use std::pow which is faster + if ( !__builtin_constant_p( base ) ) return std::pow(base, exp); + // otherwise, use constexpr-friendly implementations + else + #endif + if constexpr( std::is_integral::value ){ + return internal::PowIntImpl(base, exp); + } else return internal::PowDoubleImpl(base, exp); } // A fast (O(log p)) integer-power command. diff --git a/tests/data/DataNode.cc b/tests/data/DataNode.cc index e3d6efcaf8..7360df0ef8 100644 --- a/tests/data/DataNode.cc +++ b/tests/data/DataNode.cc @@ -254,8 +254,17 @@ TEST_CASE("Test DataStats", "[data]") { REQUIRE(data3.GetMax() == 8); REQUIRE(data3.GetVariance() == Approx(5.87654)); REQUIRE(data3.GetStandardDeviation() == Approx(2.42416)); - REQUIRE(data3.GetSkew() == Approx(-0.151045)); - REQUIRE(data3.GetKurtosis() == Approx(-1.3253830944)); + // skewness calculation: https://www.wolframalpha.com/input/?i=skewness+1%2C2%2C3%2C4%2C5%2C6%2C7%2C8%2C8 + // allow fallback to less precise constexpr Pow + const bool skew_res{ data3.GetSkew() == Approx(-0.150986) || data3.GetSkew() == Approx(-0.151045) }; + REQUIRE( skew_res ); + // kurtosis calculateion + // - https://www.wolframalpha.com/input/?i=N%5BKurtosis%5B%7B1%2C+2%2C+3%2C+4%2C+5%2C+6%2C+7%2C+8%2C+8%7D%5D%5D+-+3 + // - https://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm + // allow fallback to less precise constexpr Pow + const bool kurtosis_res{ data3.GetKurtosis() == Approx(-1.325383) || data3.GetKurtosis() == Approx(-1.3253830944) }; + REQUIRE( kurtosis_res ); + } TEST_CASE("Test histogram", "[data]") {