Question

Je recherche une solution efficace (éventuellement standard, élégante et facile à mettre en œuvre) pour multiplier des nombres relativement grands et stocker le résultat dans un ou plusieurs entiers:

Disons que j'ai deux entiers de 64 bits déclarés comme ceci:

uint64_t a = xxx, b = yyy; 

Lorsque je le fais a * b, comment puis-je détecter si l'opération entraîne un débordement et dans ce cas stocker le report quelque part?

Veuillez noter que je ne souhaite utiliser aucune bibliothèque à grand nombre car j'ai des contraintes sur la façon dont je stocke les numéros.

Était-ce utile?

La solution

1. Détection du débordement :

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

Modifier: correction de la division par 0 (merci Mark!)

2. Le calcul du report est très complexe. Une approche consiste à diviser les deux opérandes en demi-mots, puis à appliquer la multiplication longue à la moitié -words:

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

Pour constater qu'aucune des sommes partielles ne peut déborder, considérons le pire des cas:

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

Laisser B = 1 << 32. Nous avons ensuite

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

Je pense que cela fonctionnera. Du moins, il gère le cas test de Sjlver. En dehors de cela, il n'a pas encore été testé (et peut-être même pas compilé, car je n'ai plus de compilateur C ++ sous la main).

Autres conseils

L'idée est d'utiliser le fait suivant, ce qui est vrai pour l'opération intégrale:

a*b > c si et seulement si a > c/b

/ est la division intégrale ici.

Le pseudocode à vérifier en cas de dépassement de capacité pour les nombres positifs suit:

if (a > max_int64 / b), puis "flow " sinon & "ok &"; .

Pour gérer les zéros et les nombres négatifs, vous devez ajouter davantage de contrôles.

Le code C pour les non négatifs a et b suit:

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

Remarque:

18446744073709551615 == (1<<64)-1

Pour calculer le report, nous pouvons utiliser une approche permettant de fractionner un nombre en deux chiffres de 32 chiffres et de les multiplier comme nous le faisons sur le papier. Nous devons fractionner les chiffres pour éviter les débordements.

Le code suit:

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

Bien qu’il y ait eu plusieurs autres réponses à cette question, plusieurs d’entre elles ont un code qui n’a jamais été testé, et à ce jour, personne n’a comparé de manière adéquate les différentes options possibles.

Pour cette raison, j'ai écrit et testé plusieurs implémentations possibles (la dernière est basée sur ce code à partir d'OpenBSD, discuté sur Reddit ici ). Voici le 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);
}

Voici les résultats des tests que j'ai effectués avec différents compilateurs et systèmes (dans ce cas, tous les tests ont été réalisés sur OS X, mais les résultats doivent être similaires sur les systèmes BSD ou 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 |
+------------------+----------+----------+----------+----------+----------+

Sur la base de ces résultats, nous pouvons tirer quelques conclusions:

  • De toute évidence, l'approche par division, bien que simple et portable, est lente.
  • Aucune technique n'est clairement gagnante dans tous les cas.
  • Sur les compilateurs modernes, l'approche "utiliser un plus grand" est préférable si vous pouvez l'utiliser
  • Sur les anciens compilateurs, la méthode de la multiplication longue est préférable
  • Étonnamment, GCC 4.9.0 a des régressions de performance sur GCC 4.2.1 et GCC 4.2.1 a des régressions de performance sur GCC 3.3

Une version qui fonctionne aussi quand a == 0:

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

Si vous avez besoin non seulement de détecter un dépassement de capacité, mais également de capturer le report, il vaut mieux diviser vos chiffres en parties 32 bits. Le code est un cauchemar. ce qui suit n’est qu’un croquis:

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

Le problème ne concerne pas seulement les produits partiels, mais le fait que toutes les sommes peuvent déborder.

Si je devais le faire réellement, j'écrirais une routine de multiplication étendue dans le langage d'assemblage local. Par exemple, multipliez deux entiers de 64 bits pour obtenir un nombre de 128. résultat de bit, qui est stocké dans deux registres de 64 bits. Tout le matériel raisonnable fournit cette fonctionnalité dans une seule instruction de multiplication native & # 8212; elle n’est pas accessible uniquement à partir de & Nbsp; C.

C’est l’un des rares cas où la solution la plus élégante et la plus facile à programmer consiste à utiliser le langage assembleur. Mais ce n’est certainement pas portable: - (

Je travaille avec ce problème depuis quelques jours et je dois dire que cela m'a impressionné le nombre de fois où j'ai vu des gens dire que le meilleur moyen de savoir s'il y a eu un dépassement est de diviser le résultat, c'est-à-dire totalement inefficace et inutile. Le but de cette fonction est qu’elle doit être aussi rapide que possible.

Il existe deux options pour la détection de débordement:

1 & # 186; - Si possible, créez une variable de résultat deux fois plus grande que les multiplicateurs, par exemple:

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

Vous saurez immédiatement s'il y a eu un dépassement de capacité et le code est le plus rapide possible sans l'écrire dans le code machine. Selon le compilateur, ce code peut être amélioré dans le code machine.

2 & # 186; - Il est impossible de créer une variable de résultat deux fois plus grande que la variable multiplicateur: Ensuite, vous devez jouer avec si les conditions déterminent le meilleur chemin. Continuant avec l'exemple:

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

J'espère que ce code vous aidera à avoir un programme assez efficace et que le code est clair, sinon je vais faire quelques commentaires.

meilleures salutations.

Peut-être que le meilleur moyen de résoudre ce problème est d’avoir une fonction qui multiplie deux UInt64 et donne une paire de UInt64, une partie supérieure et une partie inférieure du résultat UInt128. Voici la solution, incluant une fonction, qui affiche le résultat en hexadécimal. Je suppose que vous préférez peut-être une solution C ++, mais j’ai une solution Swift qui montre comment gérer le problème:

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

Voici un exemple pour vérifier que la fonction fonctionne:

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

Voici une astuce pour détecter si la multiplication de deux entiers non signés déborde.

Nous observons que si nous multiplions un nombre binaire de N bits par un nombre binaire de M, le produit n'a pas plus de N + M bits.

Par exemple, si on nous demande de multiplier un nombre à trois bits par un nombre à vingt-neuf bits, nous savons que cela ne dépasse pas trente-deux 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;
}

Quelques tests: (sur un système 64 bits):

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

Les étapes de might_be_mul_oflow sont presque certainement plus lentes que le test de division, du moins sur les processeurs classiques utilisés dans les stations de travail, les serveurs et les appareils mobiles. Sur des puces sans bon support de division, cela pourrait être utile.

Il me semble qu’il existe un autre moyen de faire ce test de rejet précoce.

  1. Nous commençons avec une paire de nombres arng et brng qui sont initialisés à 0x7FFF...FFFF et à 1.

  2. Si a <= arng et b <= brng nous pouvons en conclure qu'il n'y a pas de dépassement.

  3. Sinon, nous décalons 0x3FFF...FFFF vers la droite, puis 3 vers la gauche, en ajoutant un bit à <=>, de sorte qu'ils soient <=> et <=>.

  4. Si <=> vaut zéro, terminez; sinon, répétez à 2.

La fonction ressemble maintenant à:

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

Si vous souhaitez simplement détecter un débordement, pourquoi ne pas convertir en double, effectuer la multiplication et si

| x | < 2 ^ 53, convertir en int64

| x | < 2 ^ 63, effectuez la multiplication en utilisant int64

sinon, produisez une erreur quelconque?

Cela semble fonctionner:

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;
}
Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
scroll top