Catch and compute overflow during multiplication of two large integers

Easy and fast with clang and gcc:

unsigned long long t a, b, result;
if (__builtin_umulll_overflow(a, b, &result)) {
    // overflow!!
}

This will use hardware support for overflow detection where available. By being compiler extensions it can even handle signed integer overflow (replace umul with smul), eventhough that is undefined behavior in C++.

The GNU Portability Library (Gnulib) contains a module intprops, which has macros that efficiently test whether arithmetic operations would overflow.

For example, if an overflow in multiplication would occur, INT_MULTIPLY_OVERFLOW (a, b) would yield 1.

There is a simple (and often very fast solution) which has not been mentioned yet. The solution is based on the fact that n-Bit times m-Bit multiplication does never overflow for a product width of n+m-bit or higher but overflows for all result widths smaller than n+m-1.

Because my old description might have been too difficult to read for some people, I try it again: What you need is checking the sum of leading-zeroes of both operands. It would be very easy to prove mathematically. Let x be n-Bit and y be m-Bit. z = x * y is k-Bit. Because the product can be n+m bit large at most it can overflow. Let's say. x*y is p-Bit long (without leading zeroes). The leading zeroes of the product are clz(x * y) = n+m - p. clz behaves similar to log, hence: clz(x * y) = clz(x) + clz(y) + c with c = either 1 or 0. (thank you for the c = 1 advice in the comment!) It overflows when k < p <= n+m <=> n+m - k > n+m - p = clz(x * y).

Now we can use this algorithm:

if max(clz(x * y)) = clz(x) + clz(y) +1 < (n+m - k)  --> overflow
if max(clz(x * y)) = clz(x) + clz(y) +1 == (n+m - k)  --> overflow if c = 0
else --> no overflow

How to check for overflow in the middle case? I assume, you have a multiplication instruction. Then we easily can use it to see the leading zeroes of the result, i.e.:

if clz(x * y / 2) == (n+m - k) <=> msb(x * y/2) == 1  --> overflow
else --> no overflow

You do the multiplication by treating x/2 as fixed point and y as normal integer:

msb(x * y/2) = msb(floor(x * y / 2))
floor(x * y/2) = floor(x/2) * y + (lsb(x) * floor(y/2)) = (x >> 1)*y + (x & 1)*(y >> 1)

(this result never overflows in case of clz(x)+clz(y)+1 == (n+m -k))

The trick is using builtins/intrinsics. In GCC it looks this way:

static inline int clz(int a) {
    if (a == 0) return 32; //only needed for x86 architecture
    return __builtin_clz(a);
}
/**@fn static inline _Bool chk_mul_ov(uint32_t f1, uint32_t f2)
 * @return one, if a 32-Bit-overflow occurs when unsigned-unsigned-multipliying f1 with f2 otherwise zero. */
static inline _Bool chk_mul_ov(uint32_t f1, uint32_t f2) {
    int lzsum = clz(f1) + clz(f2); //leading zero sum
    return
        lzsum < sizeof(f1)*8-1 || ( //if too small, overflow guaranteed
            lzsum == sizeof(f1)*8-1 && //if special case, do further check
            (int32_t)((f1 >> 1)*f2 + (f1 & 1)*(f2 >> 1)) < 0 //check product rightshifted by one
    );
}
...
    if (chk_mul_ov(f1, f2)) {
        //error handling
    }
...

Just an example for n = m = k = 32-Bit unsigned-unsigned-multiplication. You can generalize it to signed-unsigned- or signed-signed-multiplication. And even no multiple-bit-shift is required (because some microcontrollers implement one-bit-shifts only but sometimes support product divided by two with a single instruction like Atmega!). However, if no count-leading-zeroes instruction exists but long multiplication, this might not be better.

Other compilers have their own way of specifying intrinsics for CLZ operations. Compared to checking upper half of the multiplication the clz-method should scale better (in worst case) than using a highly optimized 128-Bit multiplication to check for 64-Bit overflow. Multiplication needs over linear overhead while count bits needs only linear overhead. This code worked out-of-the box for me when tried.

When your using e.g. 64 bits variables, implement 'number of significant bits' with nsb(var) = { 64 - clz(var); }.

clz(var) = count leading zeros in var, a builtin command for GCC and Clang, or probably available with inline assembly for your CPU.

Now use the fact that nsb(a * b) <= nsb(a) + nsb(b) to check for overflow. When smaller, it is always 1 smaller.

Ref GCC: Built-in Function: int __builtin_clz (unsigned int x) Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.

I was thinking about this today and stumbled upon this question, my thoughts led me to this result. TLDR, while I find it "elegant" in that it only uses a few lines of code (could easily be a one liner), and has some mild math that simplifies to something relatively simple conceptually, this is mostly "interesting" and I haven't tested it.

If you think of an unsigned integer as being a single digit with radix 2^n where n is the number of bits in the integer, then you can map those numbers to radians around the unit circle, e.g.

radians(x) = x * (2 * pi * rad / 2^n)

When the integer overflows, it is equivalent to wrapping around the circle. So calculating the carry is equivalent to calculating the number of times multiplication would wrap around the circle. To calculate the number of times we wrap around the circle we divide radians(x) by 2pi radians. e.g.

wrap(x) = radians(x) / (2*pi*rad)
        = (x * (2*pi*rad / 2^n)) / (2*pi*rad / 1)
        = (x * (2*pi*rad / 2^n)) * (1 / 2*pi*rad)
        = x * 1 / 2^n
        = x / 2^n

Which simplifies to

wrap(x) = x / 2^n

This makes sense. The number of times a number, for example, 15 with radix 10, wraps around is 15 / 10 = 1.5, or one and a half times. However, we can't use 2 digits here (assuming we are limited to a single 2^64 digit).

Say we have a * b, with radix R, we can calculate the carry with

Consider that: wrap(a * b) = a * wrap(b)
wrap(a * b) = (a * b) / R
a * wrap(b) = a * (b / R)
a * (b / R) = (a * b) / R

carry = floor(a * wrap(b))

Take for example a = 9 and b = 5, which are factors of 45 (i.e. 9 * 5 = 45).

wrap(5) = 5 / 10 = 0.5
a * wrap(5) = 9 * 0.5 = 4.5
carry = floor(9 * wrap(5)) = floor(4.5) = 4

Note that if the carry was 0, then we would not have had overflow, for example if a = 2, b=2.

In C/C++ (if the compiler and architecture supports it) we have to use long double.

Thus we have:

long double wrap = b / 18446744073709551616.0L; // this is b / 2^64
unsigned long carry = (unsigned long)(a * wrap); // floor(a * wrap(b))
bool overflow = carry > 0;
unsigned long c = a * b;

c here is the lower significant "digit", i.e. in base 10 9 * 9 = 81, carry = 8, and c = 1.

This was interesting to me in theory, so I thought I'd share it, but one major caveat is with the floating point precision in computers. Using long double, there may be rounding errors for some numbers when we calculate the wrap variable depending on how many significant digits your compiler/arch uses for long doubles, I believe it should be 20 more more to be sure. Another issue with this result, is that it may not perform as well as some of the other solutions simply by using floating points and division.

If you just want to detect overflow, how about converting to double, doing the multiplication and if

|x| < 2^53, convert to int64

|x| < 2^63, make the multiplication using int64

otherwise produce whatever error you want?

This seems to work:

int64_t safemult(int64_t a, int64_t b) {
  double dx;

  dx = (double)a * (double)b;

  if ( fabs(dx) < (double)9007199254740992 )
    return (int64_t)dx;

  if ( (double)INT64_MAX < fabs(dx) )
    return INT64_MAX;

  return a*b;
}