Captura y calcula el desbordamiento durante la multiplicación de dos enteros grandes

La idea es usar el siguiente hecho que es verdadero para la operación integral:

a*b > c si y solo si a > c/b

/ es división integral aquí.

El pseudocódigo para comprobar el desbordamiento de números positivos es el siguiente:

If (a > max_int64 / b) entonces "overflow" else "ok".

Para manejar ceros y números negativos debe agregar más comprobaciones.

El código C para a y b no negativo es el siguiente:

if (b > 0 && a > 18446744073709551615 / b) {
     // overflow handling
}; else {
    c = a * b;
}

Nota:

18446744073709551615 == (1<<64)-1

A calcule carry podemos usar el enfoque para dividir el número en dos dígitos de 32 y multiplicarlos a medida que lo hacemos en el papel. Tenemos que dividir los números para evitar el desbordamiento.

El código sigue:

// split input numbers into 32-bit digits
uint64_t a0 = a & ((1LL<<32)-1);
uint64_t a1 = a >> 32;
uint64_t b0 = b & ((1LL<<32)-1);
uint64_t b1 = b >> 32;


// The following 3 lines of code is to calculate the carry of d1
// (d1 - 32-bit second digit of result, and it can be calculated as d1=d11+d12),
// but to avoid overflow.
// Actually rewriting the following 2 lines:
// uint64_t d1 = (a0 * b0 >> 32) + a1 * b0 + a0 * b1;
// uint64_t c1 = d1 >> 32;
uint64_t d11 = a1 * b0 + (a0 * b0 >> 32); 
uint64_t d12 = a0 * b1;
uint64_t c1 = (d11 > 18446744073709551615 - d12) ? 1 : 0;

uint64_t d2 = a1 * b1 + c1;
uint64_t carry = d2; // needed carry stored here

Aunque ha habido varias otras respuestas a esta pregunta, yo varios de ellos tienen código que está completamente no probado, y hasta ahora nadie ha comparado adecuadamente las diferentes opciones posibles.

Por esa razón, escribí y probé varias implementaciones posibles (la última está basada en este código de OpenBSD, discutido en Reddit aquí). Aquí está el código:

/* Multiply with overflow checking, emulating clang's builtin function
 *
 *     __builtin_umull_overflow
 *
 * This code benchmarks five possible schemes for doing so.
 */

#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <limits.h>

#ifndef BOOL
    #define BOOL int
#endif

// Option 1, check for overflow a wider type
//    - Often fastest and the least code, especially on modern compilers
//    - When long is a 64-bit int, requires compiler support for 128-bits
//      ints (requires GCC >= 3.0 or Clang)

#if LONG_BIT > 32
    typedef __uint128_t long_overflow_t ;
#else
    typedef uint64_t long_overflow_t;
#endif

BOOL 
umull_overflow1(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
        long_overflow_t prod = (long_overflow_t)lhs * (long_overflow_t)rhs;
        *result = (unsigned long) prod;
        return (prod >> LONG_BIT) != 0;
}

// Option 2, perform long multiplication using a smaller type
//    - Sometimes the fastest (e.g., when mulitply on longs is a library
//      call).
//    - Performs at most three multiplies, and sometimes only performs one.
//    - Highly portable code; works no matter how many bits unsigned long is

BOOL 
umull_overflow2(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
        const unsigned long HALFSIZE_MAX = (1ul << LONG_BIT/2) - 1ul;
        unsigned long lhs_high = lhs >> LONG_BIT/2;
        unsigned long lhs_low  = lhs & HALFSIZE_MAX;
        unsigned long rhs_high = rhs >> LONG_BIT/2;
        unsigned long rhs_low  = rhs & HALFSIZE_MAX;

        unsigned long bot_bits = lhs_low * rhs_low;
        if (!(lhs_high || rhs_high)) {
            *result = bot_bits;
            return 0; 
        }
        BOOL overflowed = lhs_high && rhs_high;
        unsigned long mid_bits1 = lhs_low * rhs_high;
        unsigned long mid_bits2 = lhs_high * rhs_low;

        *result = bot_bits + ((mid_bits1+mid_bits2) << LONG_BIT/2);
        return overflowed || *result < bot_bits
            || (mid_bits1 >> LONG_BIT/2) != 0
            || (mid_bits2 >> LONG_BIT/2) != 0;
}

// Option 3, perform long multiplication using a smaller type (this code is
// very similar to option 2, but calculates overflow using a different but
// equivalent method).
//    - Sometimes the fastest (e.g., when mulitply on longs is a library
//      call; clang likes this code).
//    - Performs at most three multiplies, and sometimes only performs one.
//    - Highly portable code; works no matter how many bits unsigned long is

BOOL 
umull_overflow3(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
        const unsigned long HALFSIZE_MAX = (1ul << LONG_BIT/2) - 1ul;
        unsigned long lhs_high = lhs >> LONG_BIT/2;
        unsigned long lhs_low  = lhs & HALFSIZE_MAX;
        unsigned long rhs_high = rhs >> LONG_BIT/2;
        unsigned long rhs_low  = rhs & HALFSIZE_MAX;

        unsigned long lowbits = lhs_low * rhs_low;
        if (!(lhs_high || rhs_high)) {
            *result = lowbits;
            return 0; 
        }
        BOOL overflowed = lhs_high && rhs_high;
        unsigned long midbits1 = lhs_low * rhs_high;
        unsigned long midbits2 = lhs_high * rhs_low;
        unsigned long midbits  = midbits1 + midbits2;
        overflowed = overflowed || midbits < midbits1 || midbits > HALFSIZE_MAX;
        unsigned long product = lowbits + (midbits << LONG_BIT/2);
        overflowed = overflowed || product < lowbits;

        *result = product;
        return overflowed;
}

// Option 4, checks for overflow using division
//    - Checks for overflow using division
//    - Division is slow, especially if it is a library call

BOOL
umull_overflow4(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
        *result = lhs * rhs;
        return rhs > 0 && (SIZE_MAX / rhs) < lhs;
}

// Option 5, checks for overflow using division
//    - Checks for overflow using division
//    - Avoids division when the numbers are "small enough" to trivially
//      rule out overflow
//    - Division is slow, especially if it is a library call

BOOL
umull_overflow5(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
        const unsigned long MUL_NO_OVERFLOW = (1ul << LONG_BIT/2) - 1ul;
        *result = lhs * rhs;
        return (lhs >= MUL_NO_OVERFLOW || rhs >= MUL_NO_OVERFLOW) &&
            rhs > 0 && SIZE_MAX / rhs < lhs;
}

#ifndef umull_overflow
    #define umull_overflow2
#endif

/*
 * This benchmark code performs a multiply at all bit sizes, 
 * essentially assuming that sizes are logarithmically distributed.
 */

int main()
{
        unsigned long i, j, k;
        int count = 0;
        unsigned long mult;
        unsigned long total = 0;

        for (k = 0; k < 0x40000000 / LONG_BIT / LONG_BIT; ++k)
                for (i = 0; i != LONG_MAX; i = i*2+1)
                        for (j = 0; j != LONG_MAX; j = j*2+1) {
                                count += umull_overflow(i+k, j+k, &mult);
                                total += mult;
                        }
        printf("%d overflows (total %lu)\n", count, total);
}

Aquí están los resultados, probando con varios compiladores y sistemas que tengo (en este caso, todas las pruebas se realizaron en OS X, pero los resultados deben ser similares en sistemas BSD o Linux):

+------------------+----------+----------+----------+----------+----------+
|                  | Option 1 | Option 2 | Option 3 | Option 4 | Option 5 |
|                  |  BigInt  | LngMult1 | LngMult2 |   Div    |  OptDiv  |
+------------------+----------+----------+----------+----------+----------+
| Clang 3.5 i386   |    1.610 |    3.217 |    3.129 |    4.405 |    4.398 |
| GCC 4.9.0 i386   |    1.488 |    3.469 |    5.853 |    4.704 |    4.712 |
| GCC 4.2.1 i386   |    2.842 |    4.022 |    3.629 |    4.160 |    4.696 |
| GCC 4.2.1 PPC32  |    8.227 |    7.756 |    7.242 |   20.632 |   20.481 |
| GCC 3.3   PPC32  |    5.684 |    9.804 |   11.525 |   21.734 |   22.517 |
+------------------+----------+----------+----------+----------+----------+
| Clang 3.5 x86_64 |    1.584 |    2.472 |    2.449 |    9.246 |    7.280 |
| GCC 4.9 x86_64   |    1.414 |    2.623 |    4.327 |    9.047 |    7.538 |
| GCC 4.2.1 x86_64 |    2.143 |    2.618 |    2.750 |    9.510 |    7.389 |
| GCC 4.2.1 PPC64  |   13.178 |    8.994 |    8.567 |   37.504 |   29.851 |
+------------------+----------+----------+----------+----------+----------+

En base a estos resultados, podemos sacar algunas conclusiones:

• Claramente, el enfoque basado en la división, aunque simple y portátil, es lento.
• Ninguna técnica es un claro ganador en todos los casos.
• En los compiladores modernos, el enfoque use-a-larger-int es mejor, si puede usarlo
• En compiladores antiguos, el enfoque de multiplicación larga es best
• Sorprendentemente, GCC 4.9.0 tiene regresiones de rendimiento sobre GCC 4.2.1, y GCC 4.2.1 tiene regresiones de rendimiento sobre GCC 3.3

Una versión que también funciona cuando a = = 0:

    x = a * b;
    if (a != 0 && x / a != b) {
        // overflow handling
    }

Si necesita no solo detectar desbordamiento, sino también capturar el acarreo, es mejor dividir sus números en partes de 32 bits. El código es una pesadilla; lo que sigue es solo un boceto:

#include <stdint.h>

uint64_t mul(uint64_t a, uint64_t b) {
  uint32_t ah = a >> 32;
  uint32_t al = a;  // truncates: now a = al + 2**32 * ah
  uint32_t bh = b >> 32;
  uint32_t bl = b;  // truncates: now b = bl + 2**32 * bh
  // a * b = 2**64 * ah * bh + 2**32 * (ah * bl + bh * al) + al * bl
  uint64_t partial = (uint64_t) al * (uint64_t) bl;
  uint64_t mid1    = (uint64_t) ah * (uint64_t) bl;
  uint64_t mid2    = (uint64_t) al * (uint64_t) bh;
  uint64_t carry   = (uint64_t) ah * (uint64_t) bh;
  // add high parts of mid1 and mid2 to carry
  // add low parts of mid1 and mid2 to partial, carrying
  //    any carry bits into carry...
}

El problema no es solo los productos parciales, sino el hecho de que cualquiera de las sumas puede desbordarse.

Si tuviera que hacer esto de verdad, escribiría una rutina de multiplicación extendida en el lenguaje ensamblador local. Es decir, por ejemplo, multiplicar dos enteros de 64 bits para obtener una resultado de 128 bits, que se almacena en dos registros de 64 bits. Todo el hardware razonable proporciona esta funcionalidad en una sola instrucción nativa multiply-no solo es accesible desde C.

Este es uno de esos raros casos en los que la solución más elegante y fácil de programar es en realidad usar lenguaje ensamblador. Pero ciertamente no es portátil : - (

He estado trabajando con este problema estos días y tengo que decir que me ha impresionado el número de veces que he visto a personas diciendo que la mejor manera de saber si ha habido un desbordamiento es dividir el resultado, eso es totalmente ineficiente e innecesario. El punto para esta función es que debe ser lo más rápido posible.

Hay dos opciones para la detección de desbordamiento:

1º-Si es posible crear la variable de resultado dos veces más grande que los multiplicadores, para ejemplo:

struct INT32struct {INT16 high, low;};
typedef union
{
  struct INT32struct s;
  INT32 ll;
} INT32union;

INT16 mulFunction(INT16 a, INT16 b)
{
  INT32union result.ll = a * b; //32Bits result
  if(result.s.high > 0) 
      Overflow();
  return (result.s.low)
}

Sabrá inmediatamente si ha habido un desbordamiento, y el código es el más rápido posible sin escribirlo en código máquina. Dependiendo del compilador, este código se puede mejorar en código máquina.

2º-Es imposible crear una variable de resultado dos veces más grande que la variable multiplicadora: Entonces usted debe jugar con si las condiciones para determinar el mejor camino. Continuando con el ejemplo:

INT32 mulFunction(INT32 a, INT32 b)
{

  INT32union s_a.ll = abs(a);
  INT32union s_b.ll = abs(b); //32Bits result
  INT32union result;
  if(s_a.s.hi > 0 && s_b.s.hi > 0)
  {
      Overflow();
  }
  else if (s_a.s.hi > 0)
  {
      INT32union res1.ll = s_a.s.hi * s_b.s.lo;
      INT32union res2.ll = s_a.s.lo * s_b.s.lo;
      if (res1.hi == 0)
      {
          result.s.lo = res1.s.lo + res2.s.hi;
          if (result.s.hi == 0)
          {
            result.s.ll = result.s.lo << 16 + res2.s.lo;
            if ((a.s.hi >> 15) ^ (b.s.hi >> 15) == 1)
            {
                result.s.ll = -result.s.ll; 
            }
            return result.s.ll
          }else
          {
             Overflow();
          }
      }else
      {
          Overflow();
      }
  }else if (s_b.s.hi > 0)
{

   //Same code changing a with b

}else 
{
    return (s_a.lo * s_b.lo);
}
}

Espero que este código le ayude a tener una muy eficiente programa y espero que el código sea claro, si no voy a poner algunos comentarios.

Saludos cordiales.

Quizás la mejor manera de resolver este problema es tener una función, que multiplica dos UInt64 y resulta un par de UInt64, una parte superior y una parte inferior del resultado UInt128. Aquí está la solución, incluyendo una función, que muestra el resultado en hexadecimal. Supongo que quizás prefieras una solución C++, pero tengo una solución Swift que muestra cómo manejar el problema:

func hex128 (_ hi: UInt64, _ lo: UInt64) -> String
{
    var s: String = String(format: "%08X", hi >> 32)
                  + String(format: "%08X", hi & 0xFFFFFFFF)
                  + String(format: "%08X", lo >> 32)
                  + String(format: "%08X", lo & 0xFFFFFFFF)
    return (s)
}

func mul64to128 (_ multiplier: UInt64, _ multiplicand : UInt64)
             -> (result_hi: UInt64, result_lo: UInt64)
{
    let x: UInt64 = multiplier
    let x_lo: UInt64 = (x & 0xffffffff)
    let x_hi: UInt64 = x >> 32

    let y: UInt64 = multiplicand
    let y_lo: UInt64 = (y & 0xffffffff)
    let y_hi: UInt64 = y >> 32

    let mul_lo: UInt64 = (x_lo * y_lo)
    let mul_hi: UInt64 = (x_hi * y_lo) + (mul_lo >> 32)
    let mul_carry: UInt64 = (x_lo * y_hi) + (mul_hi & 0xffffffff)
    let result_hi: UInt64 = (x_hi * y_hi) + (mul_hi >> 32) + (mul_carry >> 32)
    let result_lo: UInt64 = (mul_carry << 32) + (mul_lo & 0xffffffff)

    return (result_hi, result_lo)
}

Aquí hay un ejemplo para verificar, que la función funciona:

var c: UInt64 = 0
var d: UInt64 = 0

(c, d) = mul64to128(0x1234567890123456, 0x9876543210987654)
// 0AD77D742CE3C72E45FD10D81D28D038 is the result of the above example
print(hex128(c, d))

(c, d) = mul64to128(0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF)
// FFFFFFFFFFFFFFFE0000000000000001 is the result of the above example
print(hex128(c, d))

Si solo desea detectar desbordamiento, ¿qué tal convertir a doble, haciendo la multiplicación y si

/ x /

/ x /

De lo contrario producir cualquier error que desee?

Esto parece funcionar:

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;
}

Aquí hay un truco para detectar si la multiplicación de dos enteros sin signo se desborda.

Hacemos la observación de que si multiplicamos un número binario de N-bit-wide con un número binario de M-bit-wide, el producto no tiene más de N + M bits.

Por ejemplo, si se nos pide multiplicar un número de tres bits por un número de veintinueve bits, sabemos que este no rebosa treinta y dos bits.

#include <stdlib.h>
#include <stdio.h>

int might_be_mul_oflow(unsigned long a, unsigned long b)
{
  if (!a || !b)
    return 0;

  a = a | (a >> 1) | (a >> 2) | (a >> 4) | (a >> 8) | (a >> 16) | (a >> 32);
  b = b | (b >> 1) | (b >> 2) | (b >> 4) | (b >> 8) | (b >> 16) | (b >> 32);

  for (;;) {
    unsigned long na = a << 1;
    if (na <= a)
      break;
    a = na;
  }

  return (a & b) ? 1 : 0;
}

int main(int argc, char **argv)
{
  unsigned long a, b;
  char *endptr;

  if (argc < 3) {
    printf("supply two unsigned long integers in C form\n");
    return EXIT_FAILURE;
  }

  a = strtoul(argv[1], &endptr, 0);

  if (*endptr != 0) {
    printf("%s is garbage\n", argv[1]);
    return EXIT_FAILURE;
  }

  b = strtoul(argv[2], &endptr, 0);

  if (*endptr != 0) {
    printf("%s is garbage\n", argv[2]);
    return EXIT_FAILURE;
  }

  if (might_be_mul_oflow(a, b))
    printf("might be multiplication overflow\n");

  {
    unsigned long c = a * b;
    printf("%lu * %lu = %lu\n", a, b, c);
    if (a != 0 && c / a != b)
      printf("confirmed multiplication overflow\n");
  }

  return 0;
}

Un puñado de pruebas: (en 64 bits system):

$ ./uflow 0x3 0x3FFFFFFFFFFFFFFF
3 * 4611686018427387903 = 13835058055282163709

$ ./uflow 0x7 0x3FFFFFFFFFFFFFFF
might be multiplication overflow
7 * 4611686018427387903 = 13835058055282163705
confirmed multiplication overflow

$ ./uflow 0x4 0x3FFFFFFFFFFFFFFF
might be multiplication overflow
4 * 4611686018427387903 = 18446744073709551612

$ ./uflow 0x5 0x3FFFFFFFFFFFFFFF
might be multiplication overflow
5 * 4611686018427387903 = 4611686018427387899
confirmed multiplication overflow

Los pasos en might_be_mul_oflow son casi seguramente más lentos que simplemente hacer la prueba de división, al menos en los procesadores convencionales utilizados en estaciones de trabajo de escritorio, servidores y dispositivos móviles. En fichas sin un buen soporte de división, podría ser útil.

Se me ocurre que hay otra manera de hacer esta prueba temprana de rechazo.

1. Comenzamos con un par de números arng y brng que se inicializan a 0x7FFF...FFFF y 1.

2. Si a <= arng y b <= brng podemos concluir que no hay desbordamiento.

3. De lo contrario, desplazamos arng a la derecha, y desplazamos brng a la izquierda, agregando un bit a brng, de modo que sean 0x3FFF...FFFF y 3.

4. Si arng es cero, terminar; de lo contrario repetir en 2.

La función ahora se ve como:

int might_be_mul_oflow(unsigned long a, unsigned long b)
{
  if (!a || !b)
    return 0;

  {
    unsigned long arng = ULONG_MAX >> 1;
    unsigned long brng = 1;

    while (arng != 0) {
      if (a <= arng && b <= brng)
        return 0;
      arng >>= 1;
      brng <<= 1;
      brng |= 1;
    }

    return 1;
  }
}