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 #include #include #include int main() { namespace bu = boost::units; namespace si = boost::units::si; using dimless = bu::quantity; bu::quantity dx = 1.0f * si::meter, dy = 1.0f * si::meter; bu::quantity v = 1.0f * si::meter_per_second; bu::quantity 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 ‘’) 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 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 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 #include #include #include #include template using static_power_helper = boost::units::unit< typename boost::units::static_power< typename boost::units::get_dimension<Unit>::type, typename boost::units::static_rational::type>::type, typename boost::units::get_system<Unit>::type>; template auto ratpow(boost::units::quantity const &x) -> boost::units::quantity, Y> { using std::pow; return decltype(ratpow<N>(x))::from_value(pow(x.value(), N)); } template auto ratpow(boost::units::quantity const &x) -> boost::units::quantity, Y> { using std::pow; return decltype(ratpow(x))::from_value(pow(x.value(), (Y) N / M)); } template auto msqrt(boost::units::quantity 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 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