[C++] Comment éviter le débordement dans expr. A B C D


Answers

La solution la plus simple et la plus générale consiste à utiliser une représentation qui ne peut pas déborder, soit en utilisant une bibliothèque d'entiers longs (par exemple http://gmplib.org/ ), soit en utilisant une structure ou un tableau et en implémentant une sorte de multiplication longue ( c'est-à-dire séparer chaque nombre en deux moitiés de 32 bits et effectuer la multiplication comme ci-dessous:

(R1 + R2 * 2^32 + R3 * 2^64 + R4 * 2^96) = R = A*B = (A1 + A2 * 2^32) * (B1 + B2 * 2^32) 
R1 = (A1*B1) % 2^32
R2 = ((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) % 2^32
R3 = (((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) / 2^32 + (A1*B2) / 2^32 + (A2*B1) / 2^32 + (A2*B2) % 2^32) %2^32
R4 = ((((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) / 2^32 + (A1*B2) / 2^32 + (A2*B1) / 2^32 + (A2*B2) % 2^32) / 2^32) + (A2*B2) / 2^32

En supposant que le résultat final correspond à 64 bits, vous n'avez pas vraiment besoin de la plupart des bits de R3 et aucun de R4

Question

J'ai besoin de calculer une expression qui ressemble à: A*B - C*D , où leurs types sont: signed long long int A, B, C, D; Chaque nombre peut être vraiment grand (ne pas déborder de son type). Alors que A*B peut provoquer un débordement, l'expression A*B - C*D peut être très petite en même temps. Comment puis-je le calculer correctement?

Par exemple: MAX * MAX - (MAX - 1) * (MAX + 1) == 1 , où MAX = LLONG_MAX - n et n - un certain nombre naturel.




Vous pouvez écrire chaque nombre dans un tableau, chaque élément étant un chiffre et faire les calculs en tant que polynomials . Prenez le polynôme résultant, qui est un tableau, et calculez le résultat en multipliant chaque élément du tableau par 10 à la puissance de la position dans le tableau (la première position étant la plus grande et la dernière étant zéro).

Le nombre 123 peut être exprimé comme suit:

123 = 100 * 1 + 10 * 2 + 3

pour lequel vous venez de créer un tableau [1 2 3] .

Vous faites cela pour tous les nombres A, B, C et D, et ensuite vous les multipliez comme des polynômes. Une fois que vous avez le polynôme résultant, vous venez de reconstruire le nombre de celui-ci.




Cela devrait fonctionner (je pense):

signed long long int a = 0x7ffffffffffffffd;
signed long long int b = 0x7ffffffffffffffd;
signed long long int c = 0x7ffffffffffffffc;
signed long long int d = 0x7ffffffffffffffe;
signed long long int bd = b / d;
signed long long int bdmod = b % d;
signed long long int ca = c / a;
signed long long int camod = c % a;
signed long long int x = (bd - ca) * a * d - (camod * d - bdmod * a);

Voici ma dérivation:

x = a * b - c * d
x / (a * d) = (a * b - c * d) / (a * d)
x / (a * d) = b / d - c / a

now, the integer/mod stuff:
x / (a * d) = (b / d + ( b % d ) / d) - (c / a + ( c % a ) / a )
x / (a * d) = (b / d - c / a) - ( ( c % a ) / a - ( b % d ) / d)
x = (b / d - c / a) * a * d - ( ( c % a ) * d - ( b % d ) * a)



Choisissez K = a big number (par exemple K = A - sqrt(A) )

A*B - C*D = (A-K)*(B-K) - (C-K)*(D-K) + K*(A-C+B-D); // Avoid overflow.

Pourquoi?

(A-K)*(B-K) = A*B - K*(A+B) + K^2
(C-K)*(D-K) = C*D - K*(C+D) + K^2

=>
(A-K)*(B-K) - (C-K)*(D-K) = A*B - K*(A+B) + K^2 - {C*D - K*(C+D) + K^2}
(A-K)*(B-K) - (C-K)*(D-K) = A*B - C*D - K*(A+B) + K*(C+D) + K^2 - K^2
(A-K)*(B-K) - (C-K)*(D-K) = A*B - C*D - K*(A+B-C-D)

=>
A*B - C*D = (A-K)*(B-K) - (C-K)*(D-K) + K*(A+B-C-D)

=>
A*B - C*D = (A-K)*(B-K) - (C-K)*(D-K) + K*(A-C+B-D)

Notez que parce que A, B, C et D sont de grands nombres, AC et BD sont de petits nombres.




Si vous savez que le résultat final est représentable dans votre type entier, vous pouvez effectuer ce calcul rapidement en utilisant le code ci-dessous. Étant donné que la norme C spécifie que l'arithmétique non signée est arithmétique modulo et ne déborde pas, vous pouvez utiliser un type non signé pour effectuer le calcul.

Le code suivant suppose qu'il existe un type non signé de la même largeur et que le type signé utilise tous les motifs binaires pour représenter les valeurs (aucune représentation d'interruption, le minimum du type signé est le négatif de la moitié du module du type non signé). Si cela ne tient pas dans une implémentation C, des ajustements simples peuvent être apportés à la routine ConvertToSigned pour cela.

Ce qui suit utilise signed char et unsigned char pour démontrer le code. Pour votre implémentation, changez la définition de Signed en typedef signed long long int Signed; et la définition de Unsigned to typedef unsigned long long int Unsigned; .

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


//  Define the signed and unsigned types we wish to use.
typedef signed char   Signed;
typedef unsigned char Unsigned;

//  uHalfModulus is half the modulus of the unsigned type.
static const Unsigned uHalfModulus = UCHAR_MAX/2+1;

//  sHalfModulus is the negation of half the modulus of the unsigned type.
static const Signed   sHalfModulus = -1 - (Signed) (UCHAR_MAX/2);


/*  Map the unsigned value to the signed value that is the same modulo the
    modulus of the unsigned type.  If the input x maps to a positive value, we
    simply return x.  If it maps to a negative value, we return x minus the
    modulus of the unsigned type.

    In most C implementations, this routine could simply be "return x;".
    However, this version uses several steps to convert x to a negative value
    so that overflow is avoided.
*/
static Signed ConvertToSigned(Unsigned x)
{
    /*  If x is representable in the signed type, return it.  (In some
        implementations, 
    */
    if (x < uHalfModulus)
        return x;

    /*  Otherwise, return x minus the modulus of the unsigned type, taking
        care not to overflow the signed type.
    */
    return (Signed) (x - uHalfModulus) - sHalfModulus;
}


/*  Calculate A*B - C*D given that the result is representable as a Signed
    value.
*/
static signed char Calculate(Signed A, Signed B, Signed C, Signed D)
{
    /*  Map signed values to unsigned values.  Positive values are unaltered.
        Negative values have the modulus of the unsigned type added.  Because
        we do modulo arithmetic below, adding the modulus does not change the
        final result.
    */
    Unsigned a = A;
    Unsigned b = B;
    Unsigned c = C;
    Unsigned d = D;

    //  Calculate with modulo arithmetic.
    Unsigned t = a*b - c*d;

    //  Map the unsigned value to the corresponding signed value.
    return ConvertToSigned(t);
}


int main()
{
    //  Test every combination of inputs for signed char.
    for (int A = SCHAR_MIN; A <= SCHAR_MAX; ++A)
    for (int B = SCHAR_MIN; B <= SCHAR_MAX; ++B)
    for (int C = SCHAR_MIN; C <= SCHAR_MAX; ++C)
    for (int D = SCHAR_MIN; D <= SCHAR_MAX; ++D)
    {
        //  Use int to calculate the expected result.
        int t0 = A*B - C*D;

        //  If the result is not representable in signed char, skip this case.
        if (t0 < SCHAR_MIN || SCHAR_MAX < t0)
            continue;

        //  Calculate the result with the sample code.
        int t1 = Calculate(A, B, C, D);

        //  Test the result for errors.
        if (t0 != t1)
        {
            printf("%d*%d - %d*%d = %d, but %d was returned.\n",
                A, B, C, D, t0, t1);
            exit(EXIT_FAILURE);
        }
    }
    return 0;
}



Si le résultat tient dans un long long int alors l'expression A * BC * D est correcte car il effectue le mod arithmétique 2 ^ 64, et donnera le résultat correct. Le problème est de savoir si le résultat s'inscrit dans un long long int. Pour détecter cela, vous pouvez utiliser l'astuce suivante en utilisant des doubles:

if( abs( (double)A*B - (double)C*D ) > MAX_LLONG ) 
    Overflow
else 
    return A*B-C*D;

Le problème avec cette approche est que vous êtes limité par la précision de la mantisse des doubles (54bits?) Donc vous devez limiter les produits A * B et C * D à 63 + 54 bits (ou probablement un peu moins).




Par souci d'exhaustivité, puisque personne ne l'a mentionné, certains compilateurs (par exemple GCC) vous fournissent actuellement un entier de 128 bits de nos jours.

Ainsi, une solution facile pourrait être:

(long long)((__int128)A * B - (__int128)C * D)