Skip to content

Commit

Permalink
Fix incorrect transcedental function output for scaled dimensionless …
Browse files Browse the repository at this point in the history
…types

As reported in issue nholthaus#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 <iostream>
    #include <cmath>
    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<class ScalarUnit>
    dimensionless::scalar_t exp(const ScalarUnit x) noexcept
    {
        static_assert(traits::is_dimensionless_unit<ScalarUnit>::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::unit<std::ratio<1, 1000000>, 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>, std::ratio<0, 1> >, double, units::linear_scale>) $0 = {
      units::linear_scale<double> = (m_value = -1)
    }

    # const auto c = 1.0_um / 1.0_m;
    (lldb) p c
    (const units::dimensionless::scalar_t) $0 = {
      units::linear_scale<double> = (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::unit<std::ratio<1, 1>, 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>, std::ratio<0, 1> >, double, units::linear_scale>) $0 = {
      units::linear_scale<double> = (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 Units, typename T = UNIT_LIBDEFAULT_TYPE, template<typename> class NonLinearScale = linear_scale>
    class unit_t : public NonLinearScale<T>, units::detail::_unit_t
    {
    public:
        typedef T underlying_type;

        inline constexpr underlying_type value() const noexcept
        {
            return static_cast<underlying_type>(*this);
        }

        template<typename Ty, class = std::enable_if_t<std::is_arithmetic<Ty>::value>>
        inline constexpr Ty to() const noexcept
        {
            return static_cast<Ty>(*this);
        }

        template<class Ty, std::enable_if_t<traits::is_dimensionless_unit<Units>::value && std::is_arithmetic<Ty>::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<Ty>(units::convert<Units, unit<std::ratio<1>, units::category::scalar_unit>>((*this)()));
        }
    };

    template<typename T>
    struct linear_scale
    {
        inline constexpr T operator()() const noexcept { return m_value; }
        T m_value;
    };

    template<typename T>
    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.
  • Loading branch information
Alex Wang committed May 25, 2022
1 parent ea6d126 commit 47de749
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 8 deletions.
16 changes: 8 additions & 8 deletions include/units.h
Original file line number Diff line number Diff line change
Expand Up @@ -4468,7 +4468,7 @@ namespace units
dimensionless::scalar_t exp(const ScalarUnit x) noexcept
{
static_assert(traits::is_dimensionless_unit<ScalarUnit>::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()));
}

/**
Expand All @@ -4484,7 +4484,7 @@ namespace units
dimensionless::scalar_t log(const ScalarUnit x) noexcept
{
static_assert(traits::is_dimensionless_unit<ScalarUnit>::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()));
}

/**
Expand All @@ -4499,7 +4499,7 @@ namespace units
dimensionless::scalar_t log10(const ScalarUnit x) noexcept
{
static_assert(traits::is_dimensionless_unit<ScalarUnit>::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()));
}

/**
Expand All @@ -4519,7 +4519,7 @@ namespace units
static_assert(traits::is_dimensionless_unit<ScalarUnit>::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;
}
Expand All @@ -4535,7 +4535,7 @@ namespace units
dimensionless::scalar_t exp2(const ScalarUnit x) noexcept
{
static_assert(traits::is_dimensionless_unit<ScalarUnit>::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()));
}

/**
Expand All @@ -4550,7 +4550,7 @@ namespace units
dimensionless::scalar_t expm1(const ScalarUnit x) noexcept
{
static_assert(traits::is_dimensionless_unit<ScalarUnit>::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()));
}

/**
Expand All @@ -4566,7 +4566,7 @@ namespace units
dimensionless::scalar_t log1p(const ScalarUnit x) noexcept
{
static_assert(traits::is_dimensionless_unit<ScalarUnit>::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()));
}

/**
Expand All @@ -4581,7 +4581,7 @@ namespace units
dimensionless::scalar_t log2(const ScalarUnit x) noexcept
{
static_assert(traits::is_dimensionless_unit<ScalarUnit>::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()));
}

//----------------------------------
Expand Down
20 changes: 20 additions & 0 deletions unitTests/main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -2765,18 +2765,24 @@ 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)
Expand All @@ -2786,30 +2792,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)
Expand Down

0 comments on commit 47de749

Please sign in to comment.