Cattura e calcola l'overflow durante la moltiplicazione di due numeri interi grandi

StackOverflow https://stackoverflow.com/questions/1815367

Domanda

Sto cercando una soluzione efficiente (facoltativamente standard, elegante e facile da implementare) per moltiplicare numeri relativamente grandi e archiviare il risultato in uno o più numeri interi:

Diciamo che ho due interi a 64 bit dichiarati in questo modo: 

uint64_t a = xxx, b = yyy;

Quando lo faccio a * b, come posso rilevare se l'operazione si traduce in un overflow e in questo caso conservare il carry da qualche parte?

Nota che non voglio usare nessuna libreria di grandi numeri poiché ho dei vincoli sul modo in cui memorizzo i numeri.

È stato utile?

Soluzione

1. Rilevamento dell'overflow : 

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

Modifica: risolta la divisione per 0 (grazie Marco!)

2. Il calcolo del carry è abbastanza coinvolto. Un approccio consiste nel dividere entrambi gli operandi in mezze parole, quindi applicare moltiplicazione lunga LE PAROLE: 

uint64_t hi(uint64_t x) {
    return x >> 32;
}

uint64_t lo(uint64_t x) {
    return ((1L << 32) - 1) & x;
}

void multiply(uint64_t a, uint64_t b) {
    // actually uint32_t would do, but the casting is annoying
    uint64_t s0, s1, s2, s3; 

    uint64_t x = lo(a) * lo(b);
    s0 = lo(x);

    x = hi(a) * lo(b) + hi(x);
    s1 = lo(x);
    s2 = hi(x);

    x = s1 + lo(a) * hi(b);
    s1 = lo(x);

    x = s2 + hi(a) * hi(b) + hi(x);
    s2 = lo(x);
    s3 = hi(x);

    uint64_t result = s1 << 32 | s0;
    uint64_t carry = s3 << 32 | s2;
}

Per vedere che nessuna delle somme parziali stesse può traboccare, consideriamo il caso peggiore: 

        x = s2 + hi(a) * hi(b) + hi(x)

Lascia B = 1 << 32. Abbiamo quindi 

            x <= (B - 1) + (B - 1)(B - 1) + (B - 1)
              <= B*B - 1
               < B*B

Credo che funzionerà - almeno gestisce il caso di test di Sjlver. A parte questo, non è testato (e potrebbe non essere nemmeno compilato, poiché non ho più un compilatore C ++ a portata di mano).

Altri suggerimenti

L'idea è quella di utilizzare il seguente fatto che è vero per l'operazione integrale:

a*b > c se e solo se a > c/b

/ è la divisione integrale qui.

Segue lo pseudocodice per verificare l'overflow di numeri positivi:

if (a > max_int64 / b) quindi " overflow " else " ok " .

Per gestire zero e numeri negativi è necessario aggiungere più controlli.

Segue

codice C per a e b non negativi: 

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

Nota: 

18446744073709551615 == (1<<64)-1

Per calcolare il carry possiamo usare l'approccio per dividere il numero in due 32 cifre e moltiplicarli mentre lo facciamo sul foglio. Dobbiamo dividere i numeri per evitare il trabocco.

Segue il codice: 

// 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

Anche se ci sono state molte altre risposte a questa domanda, molti di loro hanno un codice completamente non testato, e finora nessuno ha adeguatamente confrontato le diverse opzioni possibili.

Per questo motivo, ho scritto e testato diverse possibili implementazioni (l'ultima si basa su questo codice da OpenBSD, discusso su Reddit qui ). Ecco il codice: 

/* 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);
}

Ecco i risultati, test con vari compilatori e sistemi che ho (in questo caso, tutti i test sono stati eseguiti su OS X, ma i risultati dovrebbero essere simili sui sistemi 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 |
+------------------+----------+----------+----------+----------+----------+

Sulla base di questi risultati, possiamo trarre alcune conclusioni:

  • Chiaramente, l'approccio basato sulla divisione, sebbene semplice e portatile, è lento.
  • Nessuna tecnica è un chiaro vincitore in tutti i casi.
  • Sui compilatori moderni, l'approccio use-a-large-int è il migliore, se puoi usarlo
  • Sui compilatori più vecchi, l'approccio a moltiplicazione lunga è il migliore
  • Sorprendentemente, GCC 4.9.0 ha regressioni delle prestazioni su GCC 4.2.1 e GCC 4.2.1 ha regressioni delle prestazioni su GCC 3.3

Una versione che funziona anche quando a == 0: 

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

Se non è necessario solo rilevare l'overflow ma anche acquisire il carry, è meglio suddividere i numeri in parti a 32 bit. Il codice è un incubo; quello che segue è solo uno schizzo: 

#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...
}

Il problema non riguarda solo i prodotti parziali, ma il fatto che una qualsiasi delle somme può traboccare.

Se dovessi farlo davvero, scriverei una routine di moltiplicazione estesa nella lingua dell'assembly locale. Cioè, ad esempio, moltiplicherei due numeri interi a 64 bit per ottenere un 128- risultato bit, memorizzato in due registri a 64 bit. Tutto l'hardware ragionevole fornisce questa funzionalità in un'unica istruzione di moltiplicazione nativa & # 8212; non è accessibile solo da & Nbsp; C.

Questo è uno di quei rari casi in cui la soluzione più elegante e facile da programmare è effettivamente usare il linguaggio assembly. Ma non è certamente portatile :-(

Ho lavorato con questo problema in questi giorni e devo dire che mi ha colpito il numero di volte in cui ho visto persone che dicevano che il modo migliore per sapere se c'è stato un trabocco è dividere il risultato, ecco totalmente inefficiente e inutile. Il punto di questa funzione è che deve essere il più veloce possibile.

Esistono due opzioni per il rilevamento dell'overflow:

1 & # 186; - Se possibile, creare la variabile di risultato due volte più grande dei moltiplicatori, ad esempio: 

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

Saprai immediatamente se si è verificato un overflow e il codice è il più veloce possibile senza scriverlo nel codice macchina. A seconda del compilatore questo codice può essere migliorato nel codice macchina.

2 & # 186; - È impossibile creare una variabile di risultato doppia rispetto alla variabile dei moltiplicatori: Quindi dovresti giocare con le condizioni per determinare il percorso migliore. Continuando con l'esempio: 

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

Spero che questo codice ti aiuti ad avere un programma abbastanza efficiente e spero che il codice sia chiaro, in caso contrario inserirò alcuni commenti.

cordiali saluti.

Forse il modo migliore per risolvere questo problema è avere una funzione, che moltiplica due UInt64 e genera una coppia di UInt64, una parte superiore e una parte inferiore del risultato UInt128. Ecco la soluzione, inclusa una funzione, che visualizza il risultato in esadecimale. Immagino che forse preferisci una soluzione C ++, ma ho una soluzione Swift funzionante che mostra come gestire il 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)
}

Ecco un esempio per verificare che la funzione funzioni: 

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))

Ecco un trucco per rilevare se la moltiplicazione di due numeri interi senza segno trabocca.

Osserviamo che se moltiplichiamo un numero binario N-bit-wide con un numero binario M-bit-wide, il prodotto non ha più di N + M bit.

Ad esempio, se ci viene chiesto di moltiplicare un numero a tre bit con un numero a ventinove bit, sappiamo che questo non ha un overflow di trentadue bit. 

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

Una manciata di test: (su sistema a 64 bit): 

$ ./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

I passaggi in might_be_mul_oflow sono quasi certamente più lenti rispetto al semplice test di divisione, almeno sui processori tradizionali utilizzati in workstation desktop, server e dispositivi mobili. Su chip senza un buon supporto di divisione, potrebbe essere utile.

Mi viene in mente che esiste un altro modo per eseguire questo test di rifiuto anticipato.

  1. Iniziamo con una coppia di numeri arng e brng che sono inizializzati su 0x7FFF...FFFF e 1.

  2. Se a <= arng e b <= brng possiamo concludere che non c'è overflow.

  3. Altrimenti, spostiamo 0x3FFF...FFFF a destra e spostiamo 3 a sinistra, aggiungendo un bit a <=>, in modo che siano <=> e <=>.

  4. Se <=> è zero, fine; altrimenti ripeti alle 2.

La funzione ora appare come: 

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

Se vuoi solo rilevare l'overflow, che ne dici di convertirlo in doppio, fare la moltiplicazione e se

| x | lt &; 2 ^ 53, converti in int64

| x | lt &; 2 ^ 63, effettua la moltiplicazione usando int64

altrimenti produce qualunque errore tu voglia?

Questo sembra funzionare: 

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;
}
Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top