7

Recently, I saw this question which asks how you can divide integers with ceil rounding (towards positive infinity). Unfortunately the answers either don't work for signed integers or have problems with underflows and overflows.

For instance, the accepted answer has this solution:

q = 1 + ((x - 1) / y);

When x is zero, there is an underflow to ~0 and the result is incorrect.

How can you implement ceil rounding correctly for signed and unsigned integers and how do you implement other rounding modes like floor (towards negative infinity) and outwards (away from zero)?

greybeard
  • 2,249
  • 8
  • 30
  • 66
Jan Schultke
  • 17,446
  • 6
  • 47
  • 96
  • 2
    @Hi-IloveSO except that floating point numbers can't represent all integers with full precision and that this conversion takes a lot of time as compared to just using integers. [See for yourself](https://godbolt.org/z/fPrh8q). The floating point version doesn't even include the `ceil` implementation yet, [which isn't simple](https://stackoverflow.com/questions/12279914/implement-ceil-in-c). If it was simple, it would have gotten inlined. – Jan Schultke Aug 16 '20 at 12:11
  • A `double` should in most cases be able to cover the entire range of an `int`. It is simple with floating point numbers, simply add or subtract `0.5` and cast it back to an integer. The time is negligible. – Hi - I love SO Aug 16 '20 at 12:18
  • 1
    @Hi-IloveSO it might be able to cover the range of an `int`, but certainly not the range of a `uint64_t`. And once again, floating point operations take longer. Addition takes at least 3 cycles and division in the ballpark of 30-60 cycles. In this time you could round using integers many times over. This increase by `0.5` you suggesting is useful for *round* by the way, not for *floor*, *ceil*, or *up*. – Jan Schultke Aug 16 '20 at 12:21
  • 2
    @Hi-IloveSO "simply add or subtract 0.5 and cast it back to an integer" - That doesn't work correctly for all numbers. See for example https://stackoverflow.com/q/35143242/5910058 – Jesper Juhl Aug 16 '20 at 12:28
  • @JesperJuhl Interesting! Well, when it comes to division in general, it's an expensive operation that should be avoided unless working with some more sophisticated mathematical model. I'm just stating some guidelines to avoid overengineering OP's solution to whatever OP is working on. Sounds like an XY problem to me. – Hi - I love SO Aug 16 '20 at 13:17
  • 1
    @Hi-IloveSO why would you call it overengineering? The compile output is short, each function is only around 20 lines of code and it's the only solution that actually works. – Jan Schultke Aug 16 '20 at 13:32
  • Away from zero is fine, but perhaps outward. – Eric Postpischil Aug 17 '20 at 21:07

2 Answers2

15

In C++, the / division operation rounds using truncate (towards zero) by default. We can adjust the result of division towards zero to implement other rounding modes. Note that when the division has no remainder, all rounding modes are equivalent because no rounding is necessary.

With that in mind, we can implement the different rounding modes. But before we get started, we will need a helper template for the return types so that we don't use auto return types everywhere:

#include <type_traits>

/**
 * Similar to std::common_type_t<A, B>, but if A or B are signed, the result will also be signed.
 *
 * This differs from the regular type promotion rules, where signed types are promoted to unsigned types.
 */
template <typename A, typename B>
using common_signed_t =
    std::conditional_t<std::is_unsigned_v<A> && std::is_unsigned_v<B>,
                       std::common_type_t<A, B>,
                       std::common_type_t<std::make_signed_t<A>, std::make_signed_t<B>>>;

Ceil (towards +∞)

Ceil rounding is identical to truncate rounding for negative quotients, but for positive quotients and nonzero remainders we round away from zero. This means that we increment the quotient for nonzero remainders.

Thanks to if-constexpr, we can implement everything using only a single function:

template <typename Dividend, typename Divisor>
constexpr common_signed_t<Dividend, Divisor> div_ceil(Dividend x, Divisor y)
{
    if constexpr (std::is_unsigned_v<Dividend> && std::is_unsigned_v<Divisor>) {
        // quotient is always positive
        return x / y + (x % y != 0);  // uint / uint
    }
    else if constexpr (std::is_signed_v<Dividend> && std::is_unsigned_v<Divisor>) {
        auto sy = static_cast<std::make_signed_t<Divisor>>(y);
        bool quotientPositive = x >= 0;
        return x / sy + (x % sy != 0 && quotientPositive);  // int / uint
    }
    else if constexpr (std::is_unsigned_v<Dividend> && std::is_signed_v<Divisor>) {
        auto sx = static_cast<std::make_signed_t<Dividend>>(x);
        bool quotientPositive = y >= 0;
        return sx / y + (sx % y != 0 && quotientPositive);  // uint / int
    }
    else {
        bool quotientPositive = (y >= 0) == (x >= 0);
        return x / y + (x % y != 0 && quotientPositive);  // int / int
    }
}

At first glance, the implementations for signed types seem expensive, because they use both an integer division and a modulo division. However, on modern architectures division typically sets a flag that indicates whether there was a remainder, so x % y != 0 is completely free in this case.

You might also be wondering why we don't compute the quotient first and then check if the quotient is positive. This wouldn't work because we already lost precision during this division, so we can't perform this test afterwards. For example:

-1 / 2 = -0.5
// C++ already rounds towards zero
-0.5 -> 0
// Now we think that the quotient is positive, even though it is negative.
// So we mistakenly round up again:
0 -> 1

Floor (towards -∞)

Floor rounding is identical to truncate for positive quotients, but for negative quotients we round away from zero. This means that we decrement the quotient for nonzero remainders.

template <typename Dividend, typename Divisor>
constexpr common_signed_t<Dividend, Divisor> div_floor(Dividend x, Divisor y)
{
    if constexpr (std::is_unsigned_v<Dividend> && std::is_unsigned_v<Divisor>) {
        // quotient is never negative
        return x / y;  // uint / uint
    }
    else if constexpr (std::is_signed_v<Dividend> && std::is_unsigned_v<Divisor>) {
        auto sy = static_cast<std::make_signed_t<Divisor>>(y);
        bool quotientNegative = x < 0;
        return x / sy - (x % sy != 0 && quotientNegative);  // int / uint
    }
    else if constexpr (std::is_unsigned_v<Dividend> && std::is_signed_v<Divisor>) {
        auto sx = static_cast<std::make_signed_t<Dividend>>(x);
        bool quotientNegative = y < 0;
        return sx / y - (sx % y != 0 && quotientNegative);  // uint / int
    }
    else {
        bool quotientNegative = (y < 0) != (x < 0);
        return x / y - (x % y != 0 && quotientNegative);  // int / int
    }
}

The implementation is almost identical to that of div_ceil.

Away From Zero

Away from zero is the exact opposite of truncate. Basically, we need to increment or decrement depending on the sign of the quotient, but only if there is a remainder. This can be expressed as adding the sgn of the quotient onto the result:

template <typename Int>
constexpr signed char sgn(Int n)
{
    return (n > Int{0}) - (n < Int{0});
};

Using this helper function, we can fully implement up rounding:

template <typename Dividend, typename Divisor>
constexpr common_signed_t<Dividend, Divisor> div_up(Dividend x, Divisor y)
{
    if constexpr (std::is_unsigned_v<Dividend> && std::is_unsigned_v<Divisor>) {
        // sgn is always 1
        return x / y + (x % y != 0);  // uint / uint
    }
    else if constexpr (std::is_signed_v<Dividend> && std::is_unsigned_v<Divisor>) {
        auto sy = static_cast<std::make_signed_t<Divisor>>(y);
        signed char quotientSgn = sgn(x);
        return x / sy + (x % sy != 0) * quotientSgn;  // int / uint
    }
    else if constexpr (std::is_unsigned_v<Dividend> && std::is_signed_v<Divisor>) {
        auto sx = static_cast<std::make_signed_t<Dividend>>(x);
        signed char quotientSgn = sgn(y);
        return sx / y + (sx % y != 0) * quotientSgn;  // uint / int
    }
    else {
        signed char quotientSgn = sgn(x) * sgn(y);
        return x / y + (x % y != 0) * quotientSgn;  // int / int
    }
}

Unresolved Problems

Unfortunately these functions won't work for all possible inputs, which is a problem that we can not solve. For example, dividing uint32_t{3 billion} / int32_t{1} results in int32_t(3 billion) which isn't representable using a 32-bit signed integer. We get an underflow in this case.

Using larger return types would be an option for everything but 64-bit integers, where there isn't a larger alternative available. Hence, it is the responsibility of the user to ensure that when they pass an unsigned number into this function, it is equivalent to its signed representation.

Jan Schultke
  • 17,446
  • 6
  • 47
  • 96
  • 1
    `div_up()` is a confusing name, it suggest "towards +infinity"; perhaps `div_out()` or `div_away_from_0()`. 2. Can you make a single header somewhere out of all of these? e.g. pastebin, github gist, etc? – einpoklum Dec 13 '21 at 09:14
0

I would simplify and use homogeneous argument type and let the users to make explicit type casts for heterogeneous inputs, if necessary. This way possible under- and overflows are moved outside these functions. Of course normal UB cases apply, eg. divide by zero and std::numeric_limits<T>::min() divided by -1 for signed T.

#include <type_traits>

//Division round up, aka take the ceiling, aka round toward positive infinity, eg. -1.5 -> -1, 1.5 -> 2
template<typename T>
    requires std::is_integral_v<T>
constexpr T divRndUp(T a, T b) noexcept
{
    if constexpr (std::is_unsigned_v<T>)
        return a / b + (a % b != 0);
    else
        return a / b + (a % b != 0 && ((a < 0) == (b < 0)));
}

//Division round down, aka take the floor, aka round toward negative infinity, eg. -1.5 -> -2, 1.5 -> 1
template<typename T>
    requires std::is_integral_v<T>
constexpr T divRndDwn(T a, T b) noexcept
{
    if constexpr (std::is_unsigned_v<T>)
        return a / b;
    else
        return a / b - (a % b != 0 && ((a < 0) != (b < 0)));
}

//Division round out, aka round out away from zero, aka round toward infinity, eg. -1.5 -> -2, 1.5 -> 2
template<typename T>
    requires std::is_integral_v<T>
constexpr T divRndOut(T a, T b) noexcept
{
    if constexpr (std::is_unsigned_v<T>)
        return a / b + (a % b != 0);
    else
        return a / b + (a % b != 0 && ((a < 0) == (b < 0))) - (a % b != 0 && ((a < 0) != (b < 0)));
}

//Division round in, aka truncate, aka round in toward zero, aka round away from infinity, eg. -1.5 -> -1, 1.5 -> 1
template<typename T>
    requires std::is_integral_v<T>
constexpr T divRndIn(T a, T b) noexcept
{
    return a / b;
}