Überlauf während der Multiplikation zweier großer Ganzzahlen abfangen und berechnen

Ich suche nach einer effizienten (optional standardmäßigen, eleganten und einfach zu implementierenden) Lösung, um relativ große Zahlen zu multiplizieren und das Ergebnis in einer oder mehreren ganzen Zahlen zu speichern:

Nehmen wir an, ich habe zwei 64-Bit-Ganzzahlen, die wie folgt deklariert sind:

uint64_t a = xxx, b = yyy; 

Wenn ich es tue a * bwie kann ich feststellen, ob die Operation zu einem Überlauf führt, und in diesem Fall den Übertrag irgendwo speichern?

Bitte beachte, dass Ich möchte keine umfangreiche Bibliothek verwenden da ich Einschränkungen bei der Art und Weise habe, wie ich die Nummern speichere.

1. Erkennen des Überlaufs:

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

Edit: Division durch korrigiert 0 (Danke Markus!)

2. Berechnen des Übertrags ist ziemlich involviert. Ein Ansatz besteht darin, beide Operanden in Halbwörter aufzuteilen und dann anzuwenden lange Multiplikation zu den Halbwörtern:

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

uint64_t lo(uint64_t x) {
    return ((1ULL << 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;
}

Damit keine der Teilsummen selbst überlaufen kann, betrachten wir den ungünstigsten Fall:

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

Lassen B = 1 << 32. Wir haben dann

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

Ich glaube, das wird funktionieren – zumindest behandelt es Sjlvers Testfall. Abgesehen davon ist es ungetestet (und wird möglicherweise nicht einmal kompiliert, da ich keinen C++-Compiler mehr zur Hand habe).

Die Idee ist, die folgende Tatsache zu verwenden, die für den integralen Betrieb gilt:

a*b > c dann und nur dann, wenn a > c/b

/ ist hier ganzzahlige Division.

Der Pseudocode zum Prüfen auf Überlauf für positive Zahlen folgt:

if (a > max_int64 / b) then “overflow” else “ok”.

Um Nullen und negative Zahlen zu behandeln, sollten Sie weitere Prüfungen hinzufügen.

C-Code für nicht negativ a und b folgt:

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

Notiz:

18446744073709551615 == (1<<64)-1

Um den Übertrag zu berechnen, können wir den Ansatz verwenden, um die Zahl in zwei 32-stellige Zahlen aufzuteilen und sie zu multiplizieren, während wir dies auf dem Papier tun. Wir müssen Zahlen aufteilen, um einen Überlauf zu vermeiden.

Code folgt:

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

Obwohl es mehrere andere Antworten auf diese Frage gegeben hat, haben einige von ihnen Code, der völlig ungetestet ist, und bisher hat niemand die verschiedenen möglichen Optionen angemessen verglichen.

Aus diesem Grund habe ich mehrere mögliche Implementierungen geschrieben und getestet (die letzte basiert auf dieser Code von OpenBSD, diskutiert auf Reddit hier). Hier ist der Code:

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

Hier sind die Ergebnisse, Tests mit verschiedenen Compilern und Systemen, die ich habe (in diesem Fall wurden alle Tests unter OS X durchgeführt, aber die Ergebnisse sollten auf BSD- oder Linux-Systemen ähnlich sein):

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

Basierend auf diesen Ergebnissen können wir einige Schlussfolgerungen ziehen:

• Der divisionsbasierte Ansatz ist eindeutig langsam, obwohl er einfach und übertragbar ist.
• Keine Technik ist in allen Fällen ein klarer Gewinner.
• Auf modernen Compilern ist der use-a-larger-int-Ansatz am besten, wenn Sie ihn verwenden können
• Bei älteren Compilern ist der Long-Multiplication-Ansatz am besten
• Überraschenderweise weist GCC 4.9.0 Leistungsrückgänge gegenüber GCC 4.2.1 und GCC 4.2.1 Leistungsrückgänge gegenüber GCC 3.3 auf

Eine Version, die auch funktioniert, wenn a == 0:

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

Einfach und schnell mit clang und gcc:

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

Dadurch wird Hardwareunterstützung für die Überlauferkennung verwendet, sofern verfügbar. Da es sich um Compiler-Erweiterungen handelt, kann es sogar einen Überlauf von vorzeichenbehafteten Ganzzahlen verarbeiten (umul durch smul ersetzen), obwohl dies ein undefiniertes Verhalten in C++ ist.

Wenn Sie nicht nur den Überlauf erkennen, sondern auch den Übertrag erfassen müssen, zerlegen Sie Ihre Zahlen am besten in 32-Bit-Teile. Der Code ist ein Alptraum; Folgendes ist nur eine Skizze:

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

Das Problem sind nicht nur die Teilprodukte, sondern die Tatsache, dass jede der Summen überlaufen kann.

Wenn ich das wirklich tun müsste, würde ich eine Extended-Multiply-Routine in der lokalen Assemblersprache schreiben. Das heißt, multiplizieren Sie beispielsweise zwei 64-Bit-Ganzzahlen, um ein 128-Bit-Ergebnis zu erhalten, das in zwei 64-Bit-Registern gespeichert wird. Jede vernünftige Hardware bietet diese Funktionalität in einer einzigen nativen Multiplizieranweisung – sie ist nicht nur von C aus zugänglich.

Dies ist einer der seltenen Fälle, in denen die eleganteste und am einfachsten zu programmierende Lösung tatsächlich die Verwendung der Assemblersprache ist. Aber es ist sicherlich nicht tragbar 🙁

Die GNU Portability Library (Gnulib) enthält ein Modul intpropsdas über Makros verfügt, die effizient testen, ob arithmetische Operationen überlaufen würden.

Wenn zum Beispiel ein Überlauf bei der Multiplikation auftreten würde, INT_MULTIPLY_OVERFLOW (a, b) ergäbe 1.