Add support for easier math expressions in mp-units/math.h

[Becheler]

Time for 2 more suggestions! <3

I began to adapt this utility:

      template<QuantityOf<isq::length> Distance>
      struct lognormal
      {
        /// @brief Dispersal location kernel parameter
        struct param_type
        {
          /// @brief Scale parameter homogeneous to a distance
          Distance a;

          /// @brief Shape parameter determining the shape of the curve (weight of long-distance events).
          double b;
        };

        /// @brief Probability density function
        /// @param r The distance radius from the source to the target
        /// @param p Parameters of the distribution
        /// @return The value of the expression \f$ \frac{1}{ (2\pi)^{3/2}br^2} exp( - \frac{log(r/a)^2}{2b^2} ) \f$
        static constexpr auto pdf(Distance r, param_type const& p)
        {
          Distance a = p.a;
          double b = p.b;
          assert(a > 0. * m && b > 0. && r >= 0. * m);
          return (1./(std::pow(2.*M_PI, 3./.2)*b*r*r)) * mp_units::exp(- (mp_units::log(r/a)*mp_units::log(r/a) / (2.*b*b) ) );
        }
      };

But I encountered two problems:

1. I had to pinpoint the library to use: std::pow when double are manipulated (mp_units::pow does not allow dimensionless types), and to use mp_units::pow when dimensions quantities are manipulated. This mix made the code a bit weirder than it should be, as I expected to be able to use mp_units::pow everywhere.
2. The mp_unit::log function seems to have been forgotten during the implementation ;)

I have been super happy with v2.0 so far! Great great job! I don't know exactly why but I find it easier to use than v1 :p

[JohelEGP]

You need to explicitly mark b as a dimensionless quantity with q * one.

[Becheler]

Neat! Thanks! Although it will not solve the non-existing mp_units::log problem, right?
On a side note, i am under the impression that mp-units math functions are not constexpr, correct ?

[JohelEGP]

It seems like some are marked constexpr, and others not.

[mpusz]

I had to pinpoint the library to use: std::pow when double are manipulated (mp_units::pow does not allow dimensionless types), and to use mp_units::pow when dimensions quantities are manipulated. This mix made the code a bit weirder than it should be, as I expected to be able to use mp_units::pow everywhere.

You should actually use pow everywhere in an unqualified way. The correct namespace should be found with ADL. Additionally, as we do not have CPOs for math functions (see https://youtu.be/eUdz0WvOMm0?si=jGXAHLXXk97ejNJf&t=4332) you should do using std::pow; to bring the default implementation for fundamental types.

[mpusz]

The mp_unit::log function seems to have been forgotten during the implementation ;)

math.h is known to be incomplete. We add things there one by one when they are needed. Please submit the PR with log and possibly all other functions that you might need.

When we have more time, we will probably extend the math.h but for now we have other priorities to scope on.

[mpusz]

On a side note, i am under the impression that mp-units math functions are not constexpr, correct ?

All math functions aside the trigonometric ones are constexpr in math.h but they might not work at compile time for most of the compilers as the std ones are getting constexpr only from C++23 and C++26. See: https://wg21.link/P0533R9 and https://wg21.link/P1383R2.

[Becheler]

Amazing! Thanks for the explanations 👍🏽 I will try to put together a PR 🤔

[Becheler]

Hi mp-units team!
I am still having some problems regarding the math function. I tried marking b as dimensionless with *one as in :

 template <QuantityOf<isq::length> Distance = mp_units::quantity<mp_units::si::metre>> struct MyStruct {
    // ...
    static constexpr auto pdf(Distance r, param_type const &p)
    {
        Distance a = p.a;
        double b = p.b;
        assert(a > 0. * m && b > 0. && r >= 0. * m);
        return b / (2. * std::numbers::pi * a * a * tgamma(2. / b)) * exp(- pow(r, b*one) / pow(a, b*one));
    }
 };

But I get a error: no matching function for call to 'pow(mp_units::quantity<mp_units::si::metre()>&, mp_units::quantity<mp_units::one(), double>)'

[JohelEGP]

mp_units::pow takes the exponent ratio as a pair of template arguments.

[mpusz]

@Becheler, it will always help us to answer you faster and more accurately if you provide a short repro in the Compiler Explorer.

[mpusz]

Calling pow() on a quantity changes its type. This is why an exponent can't be a runtime variable. We have to know it at compile time to return a correct type (e.g., length, area, volume, etc.).

[Becheler]

@mpusz @JohelEGP Indeed with a CE example it makes more sense.
Here it is : https://godbolt.org/z/cbvsezYos

[mpusz]

Maybe you wanted to just do pow() on a numerical value of the quantity and not on the quantity itself? I do not know if it makes sense, but in such a case the quantity type will not change.

[JohelEGP]

#483 rears it head?

[mpusz]

It's not just about the unit but, first and foremost, about the quantity type.

[JohelEGP]

Yeah, that's the first bullet of #483 (comment).

[Becheler]

Thanks guys!

At least in this expression, calling pow on the numerical values of the quantites does not change the dimension of the resulting expression (1/m^2). I will have to look at others distributions to make sure of it though (what defeats a bit the automation of dimension analysis but oh well, mp_units rocks as it is! 🥇

Here is a working example: https://godbolt.org/z/5n8dreTez