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