How to use sqrt and ceil with Boost::multiprecision?

Question

Do you know how to do this simple line of code without error using Boost::multiprecison ?

boost::multiprecision::cpp_int v, uMax, candidate;
//...
v += 6 * ceil((sqrt(uMax * uMax - candidate) - v) / 6);

Using MSVC there is an error for "sqrt" and it's possible to fix it with:

v += 6 * ceil((sqrt(static_cast<boost::multiprecision::cpp_int>(uMax * uMax - candidate)) - v) / 6);

Then there is an error for "ceil" and it's possible to fix it with:

namespace bmp = boost::multiprecision;
typedef bmp::number<bmp::cpp_dec_float<0>> float_bmp;
v += 6 * ceil(static_cast<float_bmp>((sqrt(static_cast<bmp::cpp_int>(uMax * uMax - candidate)) - v) / 6));

Then there is an error of "generic interconvertion" !?!

I think there should be a more elegant way to realize a so simple line of code, isn't it? Let me know if you have some ideas about it please.

Regards.

Haematoma

Answer 1

The "problem" (it's actually a feature) is that you are using the number<> frontend with template expressions enabled.

This means that many operations can be greatly optimized or even eliminated before code is generated by the compiler.

You have two options:

1. break things down

   using BF = boost::multiprecision::cpp_bin_float_100;
   using BI = boost::multiprecision::cpp_int;
   BI v = 1, uMax = 9, candidate = 1;
   
   //v += 6 * ceil((sqrt(uMax * uMax - candidate) - v) / 6);
   BF tmp1(uMax * uMax - candidate);
   BF tmp2(sqrt(tmp1) - BF(v));
   BF tmp3(ceil(tmp2 / 6));
   BI tmp4(tmp3.convert_to<BI>());
   std::cout << tmp1 << " " << tmp2 << " " << tmp3 << " " << tmp4 << "\n";
   
   v = v + 6*tmp4;
   

   So you could write

   v += 6*ceil((sqrt(BF(uMax * uMax - candidate)) - BF(v)) / 6).convert_to<BI>();
   

   Which works by forcing evaluation of expression templates (and potentially lossy conversion from float -> integer using convert_to<>).

2. In general you could switch to non-expression-template versions of the types:

   using BF = mp::number<mp::cpp_bin_float_100::backend_type, mp::et_off>;
   using BI = mp::number<mp::cpp_int::backend_type, mp::et_off>;
   

   In this particular case it doesn't change much because you still have to do type "coercions" from integer -> float -> integer:

   v += 6 * ceil((sqrt(BF(uMax * uMax - candidate)) - BF(v)) / 6).convert_to<BI>();
   

3. By simplifying, if you make all types float instead (e.g. cpp_dec_float) you can get rid of these complicating artefacts:

   using BF = mp::number<mp::cpp_dec_float_100::backend_type, mp::et_off>;
   BF v = 1, uMax = 9, candidate = 1;
   
   v += 6 * ceil((sqrt(uMax * uMax - candidate) - v) / 6);
   

   CAVEAT Use your profiler to see that using et_off doesn't cause a performance problem on your code-base

Here's a demo program showing all three approaches:

Live On Coliru

#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/number.hpp>

int main() {
    namespace mp = boost::multiprecision;
    //v += 6 * ceil((sqrt(uMax * uMax - candidate) - v) / 6);
    {
        using BF = mp::cpp_bin_float_100;
        using BI = mp::cpp_int;
        BI v = 1, uMax = 9, candidate = 1;

#ifdef DEBUG
        BF tmp1(uMax * uMax - candidate);
        BF tmp2(sqrt(BF(uMax * uMax - candidate)) - BF(v));
        BF tmp3(ceil(tmp2 / 6));
        BI tmp4(tmp3.convert_to<BI>());
        std::cout << tmp1 << " " << tmp2 << " " << tmp3 << " " << tmp4 << "\n";
#endif

        v += 6*ceil((sqrt(BF(uMax * uMax - candidate)) - BF(v)) / 6).convert_to<BI>();
    }

    {
        using BF = mp::number<mp::cpp_bin_float_100::backend_type, mp::et_off>;
        using BI = mp::number<mp::cpp_int::backend_type, mp::et_off>;
        BI v = 1, uMax = 9, candidate = 1;

        v += 6 * ceil((sqrt(BF(uMax * uMax - candidate)) - BF(v)) / 6).convert_to<BI>();
    }

    {
        using BF = mp::number<mp::cpp_dec_float_100::backend_type, mp::et_off>;
        BF v = 1, uMax = 9, candidate = 1;

        v += 6 * ceil((sqrt(uMax * uMax - candidate) - v) / 6);
    }
}

Azoth

It's ok for the compilation errors :-) but I don't understand why ceil(BF(18)/6) is 3 and ceil(BF(18)/BF(6)) is 4. Could you help me again? – Cowles

typedef bmp::number<bmp::cpp_dec_float<0>, bmp::et_off> BF; but it's ok with typedef bmp::number<bmp::cpp_bin_float<100>, bmp::et_off> BF;. Where could I find more explanations about how to choose the best type depending on my goals? Boost 1.57 doc is not very clear for me – Cowles

My goal is to define something homogenous between integers and floats. Currently I'm using namespace bmp = boost::multiprecision; and typedef bmp::number<bmp::cpp_int_backend<>, bmp::et_off> int_bmp; // arbitrary precision integer and typedef bmp::number<bmp::cpp_bin_float<100>, bmp::et_off> float_bmp; // arbitrary precision float number but I think that the float_bmp type is not the best one... – Cowles

Just use number<> with et on. It's homogeneous. It just stops being so on (potentially) lossy conversions, because that's actually doing different things. You can always write .convert_to<> in your code to get around this. – Azoth

Regarding the different outcomes (3 vs 4) this appears to be a bug, IMO: see this comparison matrix I made which shows different backends, compilers and expression-template modes. cpp_dec_float is the only type that gets wrong results (except when ET of off and ceil(mp_t(18)/6) is evaluated, but only on GNU c++ compiler) – Azoth

Thank you very much Sehe, it's a pleasure to have so helpful replies! – Cowles

I'll file the bug unless I find the cause myself soon. I haven't tried development branch of Boost Multiprecision anyways – Azoth