From 1ec32bebf5d1176e730ec596f72f6fab48b113a9 Mon Sep 17 00:00:00 2001 From: ts826848 Date: Mon, 26 Sep 2022 13:53:29 +0000 Subject: [PATCH] Fix incorrect math function output for scaled dimensionless types (#288) * Fix incorrect transcedental function output for scaled dimensionless types As reported in issue #284, some functions produce incorrect results when passed certain dimensionless quantites. For example, this program from the issue, slightly modified for brevity and to try to make a later point clearer: #include "units.h" #include #include using namespace units::literals; int main() { const auto c = 1.0_um * (-1 / 1.0_m); std::cout << "c : " << c << std::endl ; std::cout << "units::math::exp(c) : " << units::math::exp(c) << std::endl ; std::cout << "std::exp(c.value()) : " << std::exp(c.value()) << std::endl ; return EXIT_SUCCESS ; } Outputs: c : -1e-06 units::math::exp(c) : 0.367879 std::exp(c.value()) : 0.999999 value() is basically to<>(), but with the template parameter hard-coded to underlying_type, and to<>() is one of the documented ways to pass something to a non-unit-enabled API, so getting different results this way is rather surprising. The proximate cause of this difference is the use of operator()() instead of value() or to<>() by (some?) functions in units::math. In the case of the example above: template dimensionless::scalar_t exp(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); return dimensionless::scalar_t(std::exp(x())); // operator()() instead of value()/to() ^^ } The use of operator()() means that std::exp() is fed a different input than when value() is used, as demonstrated by the output from an appropriately modified version of the example program: c : -1e-06 c() : -1 c.value() : -1e-06 units::math::exp(c) : 0.367879 std::exp(c.value()) : 0.999999 However, this problem does not for all possible inputs. If the expression for c is changed as such: // const auto c = 1.0_um * (-1 / 1.0_m); const auto c = -1.0_um / 1.0_m; The issue appears to vanish: c : -1e-06 c() : -1e-06 c.value() : -1e-06 units::math::exp(c) : 0.999999 std::exp(c.value()) : 0.999999 I believe this is because the different calculations result in objects with different types: # const auto c = 1.0_um * (-1 / 1.0_m); (lldb) p c (const units::unit_t, units::base_unit, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1> >, std::ratio<0, 1>, std::ratio<0, 1> >, double, units::linear_scale>) $0 = { units::linear_scale = (m_value = -1) } # const auto c = 1.0_um / 1.0_m; (lldb) p c (const units::dimensionless::scalar_t) $0 = { units::linear_scale = (m_value = 9.9999999999999995E-7) } Where the latter is equivalent to: # const auto c = 1.0_um / 1.0_m; (lldb) p c (const units::unit_t, units::base_unit, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1>, std::ratio<0, 1> >, std::ratio<0, 1>, std::ratio<0, 1> >, double, units::linear_scale>) $0 = { units::linear_scale = (m_value = 9.9999999999999995E-7) } The difference in the results amounts to the "new" calculation having the "correct" value directly in m_value, whereas the "old" calculation has some information in the type, which could be described as the exponent in the type and the mantissa in m_value. Floating-point errors aside, they are equal. This points to the ultimate cause of the issue: operator()() and value()/to<>() appear to consider different information when producing their results. operator()() fully discards information encoded in the type, directly returning the "raw" underlying value (modified for the decibel scale if needed), while for dimensionless types value() and to<>() account for both information encoded in the type and the "raw" underlying value in their return values. This is apparent in their implementations, where operator()() has no awareness of the units::unit tag, but value()/to<>() delegate to the conversion operator for dimensionless quantities, which explicitly "normalizes" the underlying value using units::convert<>() before returning it: template class NonLinearScale = linear_scale> class unit_t : public NonLinearScale, units::detail::_unit_t { public: typedef T underlying_type; inline constexpr underlying_type value() const noexcept { return static_cast(*this); } template::value>> inline constexpr Ty to() const noexcept { return static_cast(*this); } template::value && std::is_arithmetic::value, int> = 0> inline constexpr operator Ty() const noexcept { // this conversion also resolves any PI exponents, by converting from a non-zero PI ratio to a zero-pi ratio. return static_cast(units::convert, units::category::scalar_unit>>((*this)())); } }; template struct linear_scale { inline constexpr T operator()() const noexcept { return m_value; } T m_value; }; template struct decibel_scale { inline constexpr T operator()() const noexcept { return 10 * std::log10(m_value); } T m_value; }; As a result, the problem is clear: the functions using operator()() to obtain a value to pass to std::math functions are incorrect and should be changed to use value() so relevant information is not thrown away. The functions changed in this commit are the functions grouped under the "TRANSCEDENTAL FUNCTIONS" header in the header file. This is why modf() is changed here despite not being a transcedental function. Some other functions suffer from a similar issue, but those will be addressed in different commits. * Fix incorrect trig function output for scaled dimensionless types Similar to the previous commit, but for certain trigonometric functions. I don't think that explicitly using value() is strictly necessary for the other trig functions since they already use convert<>() to normalize, so I don't think normalization through value() is necessary. I have not tested it, though. Not entirely sure the modifications to the test cases are the best way to make the desired changes. The idea is to use the same inputs/outputs, but to change the unit type to have a conversion factor not equal to 1. Co-authored-by: Alex Wang --- include/units.h | 28 +++++++-------- unitTests/main.cpp | 86 +++++++++++++++++++++++++++++++++++++++------- 2 files changed, 88 insertions(+), 26 deletions(-) diff --git a/include/units.h b/include/units.h index 42f6f1d6..b3fd9531 100644 --- a/include/units.h +++ b/include/units.h @@ -4283,7 +4283,7 @@ namespace units angle::radian_t acos(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return angle::radian_t(std::acos(x())); + return angle::radian_t(std::acos(x.value())); } #endif @@ -4299,7 +4299,7 @@ namespace units angle::radian_t asin(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return angle::radian_t(std::asin(x())); + return angle::radian_t(std::asin(x.value())); } #endif @@ -4319,7 +4319,7 @@ namespace units angle::radian_t atan(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return angle::radian_t(std::atan(x())); + return angle::radian_t(std::atan(x.value())); } #endif @@ -4410,7 +4410,7 @@ namespace units angle::radian_t acosh(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return angle::radian_t(std::acosh(x())); + return angle::radian_t(std::acosh(x.value())); } #endif @@ -4426,7 +4426,7 @@ namespace units angle::radian_t asinh(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return angle::radian_t(std::asinh(x())); + return angle::radian_t(std::asinh(x.value())); } #endif @@ -4444,7 +4444,7 @@ namespace units angle::radian_t atanh(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return angle::radian_t(std::atanh(x())); + return angle::radian_t(std::atanh(x.value())); } #endif @@ -4468,7 +4468,7 @@ namespace units dimensionless::scalar_t exp(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return dimensionless::scalar_t(std::exp(x())); + return dimensionless::scalar_t(std::exp(x.value())); } /** @@ -4484,7 +4484,7 @@ namespace units dimensionless::scalar_t log(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return dimensionless::scalar_t(std::log(x())); + return dimensionless::scalar_t(std::log(x.value())); } /** @@ -4499,7 +4499,7 @@ namespace units dimensionless::scalar_t log10(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return dimensionless::scalar_t(std::log10(x())); + return dimensionless::scalar_t(std::log10(x.value())); } /** @@ -4519,7 +4519,7 @@ namespace units static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); UNIT_LIB_DEFAULT_TYPE intp; - dimensionless::scalar_t fracpart = dimensionless::scalar_t(std::modf(x(), &intp)); + dimensionless::scalar_t fracpart = dimensionless::scalar_t(std::modf(x.value(), &intp)); *intpart = intp; return fracpart; } @@ -4535,7 +4535,7 @@ namespace units dimensionless::scalar_t exp2(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return dimensionless::scalar_t(std::exp2(x())); + return dimensionless::scalar_t(std::exp2(x.value())); } /** @@ -4550,7 +4550,7 @@ namespace units dimensionless::scalar_t expm1(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return dimensionless::scalar_t(std::expm1(x())); + return dimensionless::scalar_t(std::expm1(x.value())); } /** @@ -4566,7 +4566,7 @@ namespace units dimensionless::scalar_t log1p(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return dimensionless::scalar_t(std::log1p(x())); + return dimensionless::scalar_t(std::log1p(x.value())); } /** @@ -4581,7 +4581,7 @@ namespace units dimensionless::scalar_t log2(const ScalarUnit x) noexcept { static_assert(traits::is_dimensionless_unit::value, "Type `ScalarUnit` must be a dimensionless unit derived from `unit_t`."); - return dimensionless::scalar_t(std::log2(x())); + return dimensionless::scalar_t(std::log2(x.value())); } //---------------------------------- diff --git a/unitTests/main.cpp b/unitTests/main.cpp index 0bbc261e..96befa34 100644 --- a/unitTests/main.cpp +++ b/unitTests/main.cpp @@ -2690,22 +2690,46 @@ TEST_F(UnitMath, tan) TEST_F(UnitMath, acos) { EXPECT_TRUE((std::is_same::type, typename std::decay::type>::value)); - EXPECT_NEAR(angle::radian_t(2).to(), acos(scalar_t(-0.41614683654)).to(), 5.0e-11); - EXPECT_NEAR(angle::degree_t(135).to(), angle::degree_t(acos(scalar_t(-0.70710678118654752440084436210485))).to(), 5.0e-12); + auto in1 = -0.41614683654; + auto in2 = -0.70710678118654752440084436210485; + auto out1 = 2; + auto out2 = 135; + EXPECT_NEAR(angle::radian_t(out1).to(), acos(scalar_t(in1)).to(), 5.0e-11); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(acos(scalar_t(in2))).to(), 5.0e-12); + auto uin1 = in1 * 1.0_m * (1.0 / (1000.0_mm)); + auto uin2 = in2 * 1.0_m * (1.0 / (1000.0_mm)); + EXPECT_NEAR(angle::radian_t(out1).to(), acos(uin1).to(), 5.0e-11); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(acos(uin2)).to(), 5.0e-12); } TEST_F(UnitMath, asin) { EXPECT_TRUE((std::is_same::type, typename std::decay::type>::value)); - EXPECT_NEAR(angle::radian_t(1.14159265).to(), asin(scalar_t(0.90929742682)).to(), 5.0e-9); - EXPECT_NEAR(angle::degree_t(45).to(), angle::degree_t(asin(scalar_t(0.70710678118654752440084436210485))).to(), 5.0e-12); + auto in1 = 0.90929742682; + auto in2 = 0.70710678118654752440084436210485; + auto out1 = 1.14159265; + auto out2 = 45; + EXPECT_NEAR(angle::radian_t(out1).to(), asin(scalar_t(in1)).to(), 5.0e-9); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(asin(scalar_t(in2))).to(), 5.0e-12); + auto uin1 = in1 * 1.0_m * (1.0 / (1000.0_mm)); + auto uin2 = in2 * 1.0_m * (1.0 / (1000.0_mm)); + EXPECT_NEAR(angle::radian_t(out1).to(), asin(uin1).to(), 5.0e-9); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(asin(uin2)).to(), 5.0e-12); } TEST_F(UnitMath, atan) { EXPECT_TRUE((std::is_same::type, typename std::decay::type>::value)); - EXPECT_NEAR(angle::radian_t(-1.14159265).to(), atan(scalar_t(-2.18503986326)).to(), 5.0e-9); - EXPECT_NEAR(angle::degree_t(-45).to(), angle::degree_t(atan(scalar_t(-1.0))).to(), 5.0e-12); + auto in1 = -2.18503986326; + auto in2 = -1.0; + auto out1 = -1.14159265; + auto out2 = -45; + EXPECT_NEAR(angle::radian_t(out1).to(), atan(scalar_t(in1)).to(), 5.0e-9); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(atan(scalar_t(in2))).to(), 5.0e-12); + auto uin1 = in1 * 1.0_m * (1.0 / (1000.0_mm)); + auto uin2 = in2 * 1.0_m * (1.0 / (1000.0_mm)); + EXPECT_NEAR(angle::radian_t(out1).to(), atan(uin1).to(), 5.0e-9); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(atan(uin2)).to(), 5.0e-12); } TEST_F(UnitMath, atan2) @@ -2743,40 +2767,64 @@ TEST_F(UnitMath, tanh) TEST_F(UnitMath, acosh) { EXPECT_TRUE((std::is_same::type, typename std::decay::type>::value)); - EXPECT_NEAR(angle::radian_t(1.316957896924817).to(), acosh(scalar_t(2.0)).to(), 5.0e-11); - EXPECT_NEAR(angle::degree_t(75.456129290216893).to(), angle::degree_t(acosh(scalar_t(2.0))).to(), 5.0e-12); + auto ins = 2; + auto out1 = 1.316957896924817; + auto out2 = 75.456129290216893; + EXPECT_NEAR(angle::radian_t(out1).to(), acosh(scalar_t(ins)).to(), 5.0e-11); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(acosh(scalar_t(ins))).to(), 5.0e-12); + auto uins = ins * 1.0_m * (1.0 / (1000.0_mm)); + EXPECT_NEAR(angle::radian_t(out1).to(), acosh(uins).to(), 5.0e-11); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(acosh(uins)).to(), 5.0e-12); } TEST_F(UnitMath, asinh) { EXPECT_TRUE((std::is_same::type, typename std::decay::type>::value)); - EXPECT_NEAR(angle::radian_t(1.443635475178810).to(), asinh(scalar_t(2)).to(), 5.0e-9); - EXPECT_NEAR(angle::degree_t(82.714219883108939).to(), angle::degree_t(asinh(scalar_t(2))).to(), 5.0e-12); + auto ins = 2; + auto out1 = 1.443635475178810; + auto out2 = 82.714219883108939; + EXPECT_NEAR(angle::radian_t(out1).to(), asinh(scalar_t(ins)).to(), 5.0e-9); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(asinh(scalar_t(ins))).to(), 5.0e-12); + auto uins = ins * 1.0_m * (1.0 / (1000.0_mm)); + EXPECT_NEAR(angle::radian_t(out1).to(), asinh(uins).to(), 5.0e-9); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(asinh(uins)).to(), 5.0e-12); } TEST_F(UnitMath, atanh) { EXPECT_TRUE((std::is_same::type, typename std::decay::type>::value)); - EXPECT_NEAR(angle::radian_t(0.549306144334055).to(), atanh(scalar_t(0.5)).to(), 5.0e-9); - EXPECT_NEAR(angle::degree_t(31.472923730945389).to(), angle::degree_t(atanh(scalar_t(0.5))).to(), 5.0e-12); + auto ins = 0.5; + auto out1 = 0.549306144334055; + auto out2 = 31.472923730945389; + EXPECT_NEAR(angle::radian_t(out1).to(), atanh(scalar_t(ins)).to(), 5.0e-9); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(atanh(scalar_t(ins))).to(), 5.0e-12); + auto uins = ins * 1.0_m * (1.0 / (1000.0_mm)); + EXPECT_NEAR(angle::radian_t(out1).to(), atanh(uins).to(), 5.0e-9); + EXPECT_NEAR(angle::degree_t(out2).to(), angle::degree_t(atanh(uins)).to(), 5.0e-12); } TEST_F(UnitMath, exp) { double val = 10.0; EXPECT_EQ(std::exp(val), exp(scalar_t(val))); + auto uval = 5.0_m * (2 / 1000_mm); + EXPECT_EQ(std::exp(uval.value()), units::math::exp(uval)); } TEST_F(UnitMath, log) { double val = 100.0; EXPECT_EQ(std::log(val), log(scalar_t(val))); + auto uval = 50.0_m * (2 / 1000_mm); + EXPECT_EQ(std::log(uval.value()), units::math::log(uval)); } TEST_F(UnitMath, log10) { double val = 100.0; EXPECT_EQ(std::log10(val), log10(scalar_t(val))); + auto uval = 50.0_m * (2 / 1000_mm); + EXPECT_EQ(std::log10(uval.value()), units::math::log10(uval)); } TEST_F(UnitMath, modf) @@ -2786,30 +2834,44 @@ TEST_F(UnitMath, modf) scalar_t modfr2; EXPECT_EQ(std::modf(val, &modfr1), modf(scalar_t(val), &modfr2)); EXPECT_EQ(modfr1, modfr2); + + auto uval = 50.0_m * (2 / 1000_mm); + double umodfr1; + decltype(uval) umodfr2; + EXPECT_EQ(std::modf(uval.value(), &umodfr1), units::math::modf(uval, &umodfr2)); + EXPECT_EQ(modfr1, modfr2); } TEST_F(UnitMath, exp2) { double val = 10.0; EXPECT_EQ(std::exp2(val), exp2(scalar_t(val))); + auto uval = 5.0_m * (2 / 1000_mm); + EXPECT_EQ(std::exp2(uval.value()), units::math::exp2(uval)); } TEST_F(UnitMath, expm1) { double val = 10.0; EXPECT_EQ(std::expm1(val), expm1(scalar_t(val))); + auto uval = 5.0_m * (2 / 1000_mm); + EXPECT_EQ(std::expm1(uval.value()), units::math::expm1(uval)); } TEST_F(UnitMath, log1p) { double val = 10.0; EXPECT_EQ(std::log1p(val), log1p(scalar_t(val))); + auto uval = 5.0_m * (2 / 1000_mm); + EXPECT_EQ(std::log1p(uval.value()), units::math::log1p(uval)); } TEST_F(UnitMath, log2) { double val = 10.0; EXPECT_EQ(std::log2(val), log2(scalar_t(val))); + auto uval = 5.0_m * (2 / 1000_mm); + EXPECT_EQ(std::log2(uval.value()), units::math::log2(uval)); } TEST_F(UnitMath, pow)