Odd behavior when taking a square root

[khansson]

I'm playing around with this library, and I found some odd behavior

This is to calculate the average speed of a particle in a Maxwellian distribution.

I wrote something simple to do this

#include <type_traits>
#include <mp-units/systems/international/international.h>
#include <mp-units/systems/isq/space_and_time.h>
#include <mp-units/systems/si/unit_symbols.h>
#include <iostream>
#include "mp-units/math.h"

using namespace mp_units;

auto AverageSpeed(const QuantityOf<isq::thermodynamic_temperature> auto &Temperature, const QuantityOf<isq::mass> auto &Mass)
{
        const auto boltz =  mp_units::si::si2019::boltzmann_constant;
        constexpr double pi = std::numbers::pi;
        const auto Internal = 8.0 * boltz * Temperature / (pi * Mass);

        using namespace mp_units::si::unit_symbols;
        const QuantityOf< pow<2>(isq::velocity)> auto Internal_converted = value_cast<(m * m / (s*s))>( Internal);
        return mp_units::pow<1,2>(Internal); //Change Internal to Internal_converted in order to compile
}

int main()
{
  using namespace mp_units::si::unit_symbols;
  using namespace mp_units::international::unit_symbols;

  const auto Temp = 1000 * K;
  const auto Mass = 10 * Da;

  const auto particleSpeed = AverageSpeed(Temp, Mass);

  std::cout << particleSpeed[m / s].number() << std::endl;
}

This currently fails to compile with at magnitude.h line 310

  if (exp.den != 1) {
    std::terminate();  // Rational powers not yet supported
  }

I think its trying to take a square root of a boltzmann constant which is not allowed. The type of Internal seem measured in K * k_b / Da.

However, if I replace the return statement of my AverageSpeed function with

return mp_units::pow<1,2>(Internal_converted);

then everything compiles. Is there a good argument for why the value_cast in necessary here?

[chiphogg]

I would expect to see this if you're doing a unit conversion between two units whose ratio has a rational power.

The reason is that we always compute unit conversion factors at compile time, but we can't currently compute rational powers at compile time. What we'd need is a constexpr version of powl(). I'm pretty sure this could be done; for example, the gcem library has one. I'm just not sure we'd want to take on the extra dependency of using it directly, so maybe we could roll our own at some point.

In any case, the code you've written is perfectly valid for the future version of this library which includes this feature. 🙂 In the meantime, I thought it would be better to produce a hard error rather than return an approximate (i.e., incorrect) result.

The reason you're able to get it to work is that when you do the conversion in squared speed, rather than speed, the conversion factor doesn't have any rational powers.

[mpusz]

Nearly all <cmath> functions will be constexpr in either C++23 (https://wg21.link/P0533R9) or C++26 (https://wg21.link/P1383R2), and I think that a few compilers already implement it this way. We just need a PR with a proper magnitudes implementation 😉

I also refactored your example a bit here: https://godbolt.org/z/v6Tf8o3hT.

[chiphogg]

Can we take advantage of that yet, though, and still maintain the C++20 compatibliity? WDYT --- will we need to detect constexpr compatibility, and provide a "hand rolled" alternative for those who don't have it?

[mpusz]

and provide a "hand rolled" alternative for those who don't have it

or do std::terminate like we are doing right now.

[khansson]

The reason you're able to get it to work is that when you do the conversion in squared speed, rather than speed, the conversion factor doesn't have any rational powers.

In theory you should still have to take a square root of one to get another conversion factor. I take it there is some constexpr check for that.

I'm pretty sure this could be done; for example, the gcem library has one. I'm just not sure we'd want to take on the extra dependency of using it directly, so maybe we could roll our own at some point.

I get the feeling, but taking a sqrt of something is too common to leave as "just wait until you have a C++23 compiler." That laves any volume/area conversion to a length as impossible let alone my energy to velocity conversions.

I'll try to put together an option dependency to gcem. Is it safe to say that only that place in compute_base_power needs to be changed?

[chiphogg]

The reason you're able to get it to work is that when you do the conversion in squared speed, rather than speed, the conversion factor doesn't have any rational powers.

In theory you should still have to take a square root of one to get another conversion factor. I take it there is some constexpr check for that.

There are two things going on here: a unit conversion, and a runtime square root computation. The reason it works when you do the unit conversion between the squared quantities is that you're essentially pushing all of the square-root-ness into the runtime computation.

I'm pretty sure this could be done; for example, the gcem library has one. I'm just not sure we'd want to take on the extra dependency of using it directly, so maybe we could roll our own at some point.

I get the feeling, but taking a sqrt of something is too common to leave as "just wait until you have a C++23 compiler." That laves any volume/area conversion to a length as impossible let alone my energy to velocity conversions.

I'll try to put together an option dependency to gcem. Is it safe to say that only that place in compute_base_power needs to be changed?

Yes, I expect that is the only place that would need to be changed.

[chiphogg]

I'm not expecting to be able to work on this before March at the earliest, because any spare bandwidth I might have would go towards the standard units library papers instead. That said, I do have some good news: Au has now implemented this feature! See aurora-opensource/au#213.

The strategy I used was to write a root(x, n) function which does a binary search to find the value r such that int_pow(r, n) is as close as possible to x. This is probably slow, but I don't think that matters because we'll only ever call it at compile time. And, it's guaranteed to produce the closest representable result, by definition! (Up to any floating point inaccuracies in our int_pow function, that is.)

Even though I don't have time to work on this for mp-units for the next several months, I'm still posting this here so that folks can discuss the implementation strategy if they want. Then it'll be just a matter of finding time to implement this.