[boost] [units] Compilation fails on quantity copy-construction

Subject: [boost] [units] Compilation fails on quantity copy-construction
From: Pavel Kretov (firegurafiku_at_[hidden])
Date: 2017-05-31 19:09:08

Hi all!

I'm trying to make use of Boost.Units in a simple FDTD simulation
program, but found a strange behavior of some kind. Consider the
following demonstration example where I try to calculate time step
from spatial steps according to the CFL condition [1]:

    #include <iostream>
    #include <boost/units/physical_dimensions.hpp>
    #include <boost/units/systems/si.hpp>
    #include <boost/units/cmath.hpp>
    #include <boost/units/io.hpp>

    int main() {
        namespace bu = boost::units;
        namespace si = boost::units::si;
        using dimless = bu::quantity<si::dimensionless, float>;

        bu::quantity<si::length, float>
            dx = 1.0f * si::meter,
            dy = 1.0f * si::meter;

        bu::quantity<si::velocity, float>
            v = 1.0f * si::meter_per_second;

        bu::quantity<si::time, float>
            dt = 0.5f / (v * bu::sqrt(1.0f/(dx*dx) + 1.0f/(dy*dy)));

        std::cout << dt << std::endl;
    }

The code looks okay to me, but when I try to compile it, I get a scary
error message:

    main.cpp:20:19: error: no match for ‘operator/’ (operand types
    are ‘float’ and ‘<THE_TERROR_THAT_FLAPS_IN_THE_NIGHT>’)

         dt = 0.5f / (v * bu::sqrt(1.0f/(dx*dx) + 1.0f/(dy*dy)));
              ~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

where the second operand type looks like

    boost::units::multiply_typeof_helper<
        boost::units::quantity<
            boost::units::unit<
                [{length_base_dimension, static_rational<1l>},
                 {time_base_dimension,   static_rational<-1l, 1l>}],
                boost::units::si::system,
            float>,
        boost::units::quantity<
            boost::units::unit<
                [{length_base_dimension, static_rational<-1l, 1l>}],
                boost::units::si::system,
            double>>::type

In the above expression [{}, {}, ...] is a pseudo-code for static
head-tail list of boost::units::dim<Unit, Rat> pairs. Note the
'double' type materialized from barely nowhere: the code explicitly
defines everything as 'float', it's the 'bu::sqrt' to blame for it to
appear. Note that replacing 'bu::sqrt' with 'bu::pow<1, 2>()' does
not help.

As I didn't feel happy with the above error message, I mutated the
code a bit and made it to work with:

    bu::quantity<si::time, float> dt;
    dt = dimless(0.5f) / (v * bu::sqrt(1.0f/(dx*dx) + 1.0f/(dy*dy)));

I had to both (a) add 'dimless', and (b) separate variable declaration
and definition. Doing only (a) didn't help.

Mkay, but what if I need to have a const quantity? (I'm a big fan of
making everything 'const'!). The only one thing I could do was to roll
out my own 'sqrt' and 'pow' functions. Here they are:

    #include <boost/units/get_system.hpp>
    #include <boost/units/get_dimension.hpp>
    #include <boost/units/static_rational.hpp>
    #include <boost/units/quantity.hpp>
    #include <boost/units/dimension.hpp>

    template<typename Unit, int N, int M = 1>
    using static_power_helper = boost::units::unit<
            typename boost::units::static_power<
                    typename boost::units::get_dimension<Unit>::type,
                    typename boost::units::static_rational<N, M>::type>::type,
            typename boost::units::get_system<Unit>::type>;

    template<int N, typename Unit, typename Y>
    auto ratpow(boost::units::quantity<Unit, Y> const &x)
            -> boost::units::quantity<static_power_helper<Unit, N>, Y> {

        using std::pow;
        return decltype(ratpow<N>(x))::from_value(pow(x.value(), N));
    }

    template<int N, int M, typename Unit, typename Y>
    auto ratpow(boost::units::quantity<Unit, Y> const &x)
            -> boost::units::quantity<static_power_helper<Unit, N, M>, Y> {

        using std::pow;
        return decltype(ratpow<N, M>(x))::from_value(pow(x.value(), (Y) N / M));
    }

    template<typename Unit, typename Y>
    auto msqrt(boost::units::quantity<Unit, Y> const &x)
            -> decltype(ratpow<1, 2>(x)) {

        return ratpow<1, 2>(x);
    }

Using them I was able to compile the original form of the expression
(as seen in the very first listing):

    bu::quantity<si::time, float> const
        dt = 0.5f / (v * msqrt(1.0f/(dx*dx) + 1.0f/(dy*dy)));

The questions:
  1. Is it possible to work around the problem without custom 'sqrt'?
  2. Is the described behavior intentional, -or-
     Is there an error in Boost.Units' implementation of 'sqrt'?
  3. If the 'boost::unit::sqrt' is fine, is there something wrong with
     my implementation?

――― Pavel K.

[1]: https://en.wikipedia.org/wiki/Courant-Friedrichs-Lewy_condition