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

Je dois calculer une expression qui ressemble à: A*B - C*D , où leurs types sont: signed long long int A, B, C, D; Chaque numéro peut être très gros (ne pas déborder de son type). Bien que A*B puisse provoquer un débordement, en même temps, l’expression A*B - C*D peut être très petite. 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.

Cela semble trop sortingvial, je suppose. Mais A*B est celui qui pourrait déborder.

Vous pouvez faire ce qui suit sans perdre en précision

 A*B - C*D = A(D+E) - (A+F)D = AD + AE - AD - DF = AE - DF ^smaller quantities E & F E = B - D (hence, far smaller than B) F = C - A (hence, far smaller than C) 

Cette décomposition peut être faite plus loin .
Comme @Gian l’a souligné, il faudra peut-être faire attention lors de l’opération de soustraction si le type n’est pas signé depuis longtemps.


Par exemple, avec le cas que vous avez dans la question, il suffit d’une itération,

  MAX * MAX - (MAX - 1) * (MAX + 1) ABCD E = B - D = -1 F = C - A = -1 AE - DF = {MAX * -1} - {(MAX + 1) * -1} = -MAX + MAX + 1 = 1 

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 longue multiplication ( 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 corresponde à 64 bits, vous n’avez pas vraiment besoin de la plupart des bits de R3 et de R4.

Notez que ce n’est pas standard car il repose sur un débordement signé enveloppant. (GCC a des indicateurs de compilation qui permettent cela.)

Mais si vous faites simplement tous les calculs long long , le résultat de l’application directe de la formule:
(A * B - C * D) sera précis tant que le résultat correct correspond à une long long .


Voici une solution qui repose uniquement sur le comportement défini par l’implémentation de la conversion d’un entier non signé en entier signé. Mais on peut s’attendre à ce que cela fonctionne sur presque tous les systèmes aujourd’hui.

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

Cela transforme les entrées en unsigned long long où le comportement de débordement est garanti par la norme. Revenir à un entier signé à la fin est la partie définie par l’implémentation, mais fonctionnera sur presque tous les environnements aujourd’hui.


Si vous avez besoin d’une solution plus pédante, je pense que vous devez utiliser “l’arithmétique longue”

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) 

Vous pourriez envisager de calculer un facteur commun le plus élevé pour toutes vos valeurs, puis de les diviser par ce facteur avant d’effectuer vos opérations arithmétiques, puis de multiplier à nouveau. Cela suppose toutefois qu’un tel facteur existe (par exemple, si A , B , C et D sont relativement premiers, ils n’auront pas de facteur commun).

De même, vous pouvez envisager de travailler sur des échelles de log, mais cela va être un peu effrayant, sous réserve d’une précision numérique.

Si le résultat tient dans un long long int, l’expression A * BC * D est correcte car elle exécute le mod arithmétique 2 ^ 64 et donnera le résultat correct. Le problème est de savoir si le résultat correspond à 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*BC*D; 

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

 E = max(A,B,C,D) A1 = A -E; B1 = B -E; C1 = C -E; D1 = D -E; 

puis

 A*B - C*D = (A1+E)*(B1+E)-(C1+E)(D1+E) = (A1+B1-C1-D1)*E + A1*B1 -C1*D1 

Vous pouvez écrire chaque numéro dans un tableau, chaque élément étant un chiffre et effectuer les calculs sous forme de polynômes . 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 égale à zéro).

Le nombre 123 peut être exprimé comme:

 123 = 100 * 1 + 10 * 2 + 3 

pour lequel vous créez simplement un tableau [1 2 3] .

Vous faites cela pour tous les nombres A, B, C et D, puis vous les multipliez en tant que polynômes. Une fois que vous avez le polynôme résultant, vous ne faites que reconstruire le nombre.

Alors qu’un signed long long int ne contiendra pas A*B , deux d’entre eux le seront. Ainsi, A*B pourrait être décomposé en termes d’arborescence d’exposant différent, chacun d’entre eux correspondant à un signed long long int .

 A1=A>>32; A0=A & 0xffffffff; B1=B>>32; B0=B & 0xffffffff; AB_0=A0*B0; AB_1=A0*B1+A1*B0; AB_2=A1*B1; 

Même chose pour C*D

Suivant la ligne droite, la sous-fraction pourrait également être appliquée à chaque paire d’ AB_i et de CD_i , en utilisant un bit de retenue supplémentaire (avec précision un entier sur 1 bit) pour chacun. Donc, si nous disons E = A * BC * D, vous obtenez quelque chose comme:

 E_00=AB_0-CD_0 E_01=(AB_0 > CD_0) == (AB_0 - CD_0 < 0) ? 0 : 1 // carry bit if overflow E_10=AB_1-CD_1 ... 

Nous continuons en transférant la moitié supérieure de E_10 à E_20 (décalage de 32 et append, puis effacer la moitié supérieure de E_10 ).

Maintenant, vous pouvez vous débarrasser du bit de E_11 en l'ajoutant au bon signe (obtenu de la partie sans report) à E_20 . Si cela déclenche un débordement, le résultat ne conviendrait pas non plus.

E_10 maintenant suffisamment d'espace pour prendre la moitié supérieure de E_00 (décalage, append, effacer) et le bit de E_01 .

E_10 peut être plus grand maintenant, alors nous répétons le transfert à E_20 .

À ce stade, E_20 doit devenir nul, sinon le résultat ne correspondra pas. La moitié supérieure de E_10 est vide comme résultat du transfert.

L'étape finale consiste à transférer à nouveau la moitié inférieure de E_20 dans E_10 .

Si l’attente que E=A*B+C*D corresponde aux signed long long int , nous avons maintenant

 E_20=0 E_10=0 E_00=E 

Si vous savez que le résultat final est représentable dans votre type d’entier, vous pouvez effectuer ce calcul rapidement en utilisant le code ci-dessous. Comme le standard C spécifie que l’arithmétique non signée est une 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 même largeur et que le type signé utilise tous les modèles de bits pour représenter des valeurs (pas de représentation de piège, 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.

Les éléments suivants utilisent un caractère signed char et un caractère unsigned char pour illustrer 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  #include  #include  // 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; } 

Vous pouvez essayer de diviser l’équation en composants plus petits qui ne débordent pas.

 AB - CD = [ A(B - N) - C( D - M )] + [AN - CM] = ( AK - CJ ) + ( AN - CM) where K = B - N J = D - M 

Si les composants continuent à déborder, vous pouvez les diviser récursivement en composants plus petits, puis les recombiner.

Je n’ai peut-être pas couvert tous les cas extrêmes, et je n’ai pas rigoureusement testé cela, mais cela implémente une technique dont je me souviens avoir utilisé dans les années 80 en essayant de faire des calculs sur 32 bits sur un processeur 16 bits. Essentiellement, vous divisez les 32 bits en deux unités de 16 bits et travaillez avec eux séparément.

 public class DoubleMaths { private static class SplitLong { // High half (or integral part). private final long h; // Low half. private final long l; // Split. private static final int SPLIT = (Long.SIZE / 2); // Make from an existing pair. private SplitLong(long h, long l) { // Let l overflow into h. this.h = h + (l >> SPLIT); this.l = l % (1l < < SPLIT); } public SplitLong(long v) { h = v >> SPLIT; l = v % (1l < < SPLIT); } public long longValue() { return (h << SPLIT) + l; } public SplitLong add ( SplitLong b ) { // TODO: Check for overflow. return new SplitLong ( longValue() + b.longValue() ); } public SplitLong sub ( SplitLong b ) { // TODO: Check for overflow. return new SplitLong ( longValue() - b.longValue() ); } public SplitLong mul ( SplitLong b ) { /* * eg 10 * 15 = 150 * * Divide 10 and 15 by 5 * * 2 * 3 = 5 * * Must therefore multiply up by 5 * 5 = 25 * * 5 * 25 = 150 */ long lbl = l * bl; long hbh = h * bh; long lbh = l * bh; long hbl = h * bl; return new SplitLong ( lbh + hbl, lbl + hbh ); } @Override public String toString () { return Long.toHexString(h)+"|"+Long.toHexString(l); } } // I'll use long and int but this can apply just as easily to long-long and long. // The aim is to calculate A*B - C*D without overflow. static final long A = Long.MAX_VALUE; static final long B = Long.MAX_VALUE - 1; static final long C = Long.MAX_VALUE; static final long D = Long.MAX_VALUE - 2; public static void main(String[] args) throws InterruptedException { // First do it with BigIntegers to get what the result should be. BigInteger a = BigInteger.valueOf(A); BigInteger b = BigInteger.valueOf(B); BigInteger c = BigInteger.valueOf(C); BigInteger d = BigInteger.valueOf(D); BigInteger answer = a.multiply(b).subtract(c.multiply(d)); System.out.println("A*B - C*D = "+answer+" = "+answer.toString(16)); // Make one and test its integrity. SplitLong sla = new SplitLong(A); System.out.println("A="+Long.toHexString(A)+" ("+sla.toString()+") = "+Long.toHexString(sla.longValue())); // Start small. SplitLong sl10 = new SplitLong(10); SplitLong sl15 = new SplitLong(15); SplitLong sl150 = sl10.mul(sl15); System.out.println("10="+sl10.longValue()+"("+sl10.toString()+") * 15="+sl15.longValue()+"("+sl15.toString()+") = "+sl150.longValue() + " ("+sl150.toString()+")"); // The real thing. SplitLong slb = new SplitLong(B); SplitLong slc = new SplitLong(C); SplitLong sld = new SplitLong(D); System.out.println("B="+Long.toHexString(B)+" ("+slb.toString()+") = "+Long.toHexString(slb.longValue())); System.out.println("C="+Long.toHexString(C)+" ("+slc.toString()+") = "+Long.toHexString(slc.longValue())); System.out.println("D="+Long.toHexString(D)+" ("+sld.toString()+") = "+Long.toHexString(sld.longValue())); SplitLong sanswer = sla.mul(slb).sub(slc.mul(sld)); System.out.println("A*B - C*D = "+sanswer+" = "+sanswer.longValue()); } } 

Impressions:

 A*B - C*D = 9223372036854775807 = 7fffffffffffffff A=7fffffffffffffff (7fffffff|ffffffff) = 7fffffffffffffff 10=10(0|a) * 15=15(0|f) = 150 (0|96) B=7ffffffffffffffe (7fffffff|fffffffe) = 7ffffffffffffffe C=7fffffffffffffff (7fffffff|ffffffff) = 7fffffffffffffff D=7ffffffffffffffd (7fffffff|fffffffd) = 7ffffffffffffffd A*B - C*D = 7fffffff|ffffffff = 9223372036854775807 

qui me semble fonctionner.

Je parie que certaines des subtilités telles que la recherche de débordement de signes, etc., me manquent, mais je pense que l'essentiel est là.

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

Ainsi, une solution simple pourrait être:

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

AB-CD = (AB-CD) * AC / AC = (B/CD/A)*A*C Ni B/C ni D/A ne peuvent déborder, donc calculez (B/CD/A) premier. Comme le résultat final ne débordera pas selon votre définition, vous pouvez effectuer les multiplications restantes en toute sécurité et calculer (B/CD/A)*A*C le résultat requirejs.

Notez que si votre entrée peut être extrêmement petite , le B/C ou le D/A peut déborder. Si cela est possible, des manipulations plus complexes peuvent être requirejses en fonction de l’inspection des entrées.

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

 A*B - C*D = (AK)*(BK) - (CK)*(DK) + K*(A-C+BD); // Avoid overflow. 

Pourquoi?

 (AK)*(BK) = A*B - K*(A+B) + K^2 (CK)*(DK) = C*D - K*(C+D) + K^2 => (AK)*(BK) - (CK)*(DK) = A*B - K*(A+B) + K^2 - {C*D - K*(C+D) + K^2} (AK)*(BK) - (CK)*(DK) = A*B - C*D - K*(A+B) + K*(C+D) + K^2 - K^2 (AK)*(BK) - (CK)*(DK) = A*B - C*D - K*(A+BCD) => A*B - C*D = (AK)*(BK) - (CK)*(DK) + K*(A+BCD) => A*B - C*D = (AK)*(BK) - (CK)*(DK) + K*(A-C+BD) 

Notez que, comme A, B, C et D sont des nombres importants, AC et BD sont donc des petits nombres.