Le Blog-Notes d'Olivier

When you do a multiply & divide operation, to ensure you don’t loose precision, you would start by the multiplication before the division. But in that case, it could be that the result of the multiplication overflows. The natural solution would go to use an intermediate storage of larger size.

#include <concepts>

// a template that gives an integer for a given number of bits
template <std::size_t bit> struct UInt { using type = void; };

// the mul&div using the larger integer type
template <std::unsigned_integral T>
T mulDivUsingLargerType(const T v, const T mulBy, const T divBy)
{
    using Larger = UInt< sizeof(T) * 8 * 2 >::type;
    return static_cast<T>(
            (static_cast<Larger>(v) * mulBy + divBy/2) / divBy);
}

template <> struct UInt<8> { using type = uint8_t; };
template <> struct UInt<16> { using type = uint16_t; };
template <> struct UInt<32> { using type = uint32_t; };
template <> struct UInt<64> { using type = uint64_t; };Langage du code : C++ (cpp)

This uses concepts (C++20). For previous C++ standards, same result can be achieved with type_traits

#include <type_traits>

template <typename T>
std::enable_if<std::is_integral_v<T> && !std::is_signed_v<T>, T> mulDivUsingLargerType(const T v, const T mulBy, const T divBy)Langage du code : PHP (php)

This technique has a problem when the quantities to manipulate are already at the maximum size handled by the processor. To overcome this limitation, one may want to pass by an intermediate double:

#include <cmath>

template <std::unsigned_integral T T>
T mulDivUsingDouble(const T v, const T mulBy, const T divBy)
{
    return static_cast<T>(std::round((static_cast<double>(v) * mulBy) / divBy));
}Langage du code : PHP (php)

Might be safe but inefficient (depends on the presence of a FPU).

There is however a safe approach that does not require this usage of another type:

template <std::unsigned_integral T>
T safeMulDiv(const T v, const T mulBy, const T divBy)
{
    T result = (v / divBy) * mulBy;
    result += ((v % divBy) * mulBy + divBy/2) / divBy;
    return result;
}Langage du code : PHP (php)

Made a little test program to exercise those 3 techniques and check that they bring the same results:

template <std::unsigned_integral T>
void tryNoOverflowMulDiv(const T mulBy, const T divBy)
{
    std::chrono::steady_clock::duration safeTime = {};
    std::chrono::steady_clock::duration byLargerTime = {};
    std::chrono::steady_clock::duration byDoubleTime = {};
    T i = 0;
    for (unsigned n = 0; n < 1000000; ++n)
    {
        auto t0 = std::chrono::steady_clock::now();
        auto tested = safeMulDiv(i, mulBy, divBy);
        auto t1 = std::chrono::steady_clock::now();
        auto ref = mulDivUsingLargerType(i, mulBy, divBy);
        auto t2 = std::chrono::steady_clock::now();
        auto ultimate = mulDivUsingDouble(i, mulBy, divBy);
        auto t3 = std::chrono::steady_clock::now();

        if (tested != ref || tested != ultimate) {
            std::cout << std::hex << "i=" << i << " ref=" << ref << " ultimate=" << ultimate << " tested=" << tested << std::dec << std::endl;
        }

        safeTime += t1-t0;
        byLargerTime += t2-t1;
        byDoubleTime += t3-t2;

        i = i * 3 + 1;
    }

    std::cout << "safe=" << safeTime.count() << ", byLarger=" << byLargerTime.count() << ", byDouble=" << byDoubleTime.count() << std::endl;
}
Langage du code : PHP (php)

It takes quite a while to execute for 32b integers. The results are:

NB: for uint64, using a larger type is not possible on a 64b machine. Using a double works only if the input value can fit in the double.