如何避免在expr溢出。 A B C D
我需要计算一个expression式，如下所示： A*B - C*D ，其types是： signed long long int A, B, C, D; 每个数字都可能很大（不会溢出）。 虽然A*B可能导致溢出，但同时expression式A*B - C*D可能非常小。 我如何正确计算它？
例如： MAX * MAX - (MAX - 1) * (MAX + 1) == 1 ，其中MAX = LLONG_MAX - n和n – 一些自然数。
我猜这似乎太重要了。 但A*B是可能溢出的那个。
你可以做到以下几点，而不会失去精度
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)
这个分解可以进一步完成 。
正如@Gian指出的那样，如果types是无符号long long，那么在减法操作中可能需要小心。
例如，对于您在问题中遇到的情况，只需要一次迭代，
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
最简单和最一般的解决scheme是使用一个不能溢出的表示法，或者使用一个长整型库（例如http://gmplib.org/ ），或者使用一个结构或数组来表示并实现一种长乘法（即将每个数字分成两个32位的一半，并执行如下的乘法：
(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
假设最终结果符合64位，实际上并不需要R3的大部分位，也不需要R4
请注意，这不是标准，因为它依赖于环绕签名溢出。 （GCC有编译器标志，使这个。）
但是，如果你只是long long完成所有的计算，直接应用公式的结果：
(A * B - C * D)将是正确的，只要正确的结果适合很long long 。
这是一种解决方法，它只依赖于将无符号整型转换为有符号整型的实现定义的行为。 但是，今天几乎所有的系统都可以运行。
(long long)((unsigned long long)A * B - (unsigned long long)C * D)
这将input转换为unsigned long long ，其中溢出行为被保证由标准包装。 最后回到一个有符号的整数是实现定义的部分，但是今天几乎所有的环境都可以工作。
如果你需要更多的迂腐解决scheme，我认为你必须使用“长算术”
这应该工作（我认为）：
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);
这是我的推导：
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)
你可以考虑计算所有值的最大公因子，然后在进行算术运算之前将它们除以该因子，然后再乘以。 这假定存在这样一个因素（例如，如果A ， B ， C和D碰巧是相对的，他们就不会有一个共同的因素）。
同样，你可以考虑在对数尺度上工作，但是这将会有点吓人，并且会受到数值精度的限制。
如果结果符合long long int，则expression式A * BC * D是可以的，因为它执行2 ^ 64的算术模式，并且会给出正确的结果。 问题是要知道结果是否适合很长的整数。 要检测这个，你可以使用下面的技巧使用双打：
if( abs( (double)A*B - (double)C*D ) > MAX_LLONG ) Overflow else return A*BC*D;
这种方法的问题在于你受限于双精度（54bits？）的尾数的精度，所以你需要将产品A * B和C * D限制在63 + 54位（或者可能less一些）。
E = max(A,B,C,D) A1 = A -E; B1 = B -E; C1 = C -E; D1 = D -E;
然后
A*B - C*D = (A1+E)*(B1+E)-(C1+E)(D1+E) = (A1+B1-C1-D1)*E + A1*B1 -C1*D1
您可以将每个数字写入一个数组中，每个元素都是一个数字，并将其计算为多项式 。 将得到的多项式作为一个数组，然后通过将数组中的每个元素与数组中的位置乘以10来计算结果（第一个位置最大，最后一个为零）。
数字123可以表示为：
123 = 100 * 1 + 10 * 2 + 3
为此只需创build一个数组[1 2 3] 。
你为所有的数字A，B，C和D做这个，然后你把它们乘以多项式。 得到多项式后，只需从中重新构造数字。
虽然signed long long int不会持有A*B ，但其中两个将会。 所以A*B可以分解为不同指数的树项，它们中的任何一个都适合一个带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;
相同的C*D
通过直接的方式，可以对每一对AB_i和CD_i同样使用一个额外的进位位（准确地说是一个1位的整数）来完成每个子对。 所以如果我们说E = A * BC * D，你会得到如下的东西：
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 ...
我们继续将E_10的上半部分E_10到E_20 （移位32并添加，然后擦除E_10上半部分）。
现在可以通过将正确的符号（从非进位部分获得）添加到E_20来消除进位位E_20 。 如果这触发溢出，结果也不适合。
E_10现在有足够的空间来取E_00的上半E_00 （移位，加， E_01 ）和进位位E_01 。
E_10现在可能再次变大，所以我们重复转移到E_20 。
此时， E_20必须变为零，否则结果不合适。 E_10的上半部分由于传输也是空的。
最后一步是将E_20的下半部分再次转换成E_10 。
如果期望E=A*B+C*D将符合signed long long int hold，我们现在有了
E_20=0 E_10=0 E_00=E
如果您知道最终结果可用整数types表示，则可以使用下面的代码快速执行此计算。 由于C标准规定无符号算术是模算术，不会溢出，因此可以使用无符号types来执行计算。
下面的代码假定有一个宽度相同的无符号types，并且带符号的types使用所有位模式来表示值（没有陷阱表示，有符号types的最小值是无符号types模的一半的负数）。 如果这不适用于C实现，则可以对ConvertToSigned例程进行简单的调整。
以下使用signed char和unsigned char来演示代码。 对于你的实现，把Signed的定义改为typedef signed long long int Signed; 并定义Unsigned为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; }
你可以尝试把方程分解成不会溢出的小部分。
AB - CD = [ A(B - N) - C( D - M )] + [AN - CM] = ( AK - CJ ) + ( AN - CM) where K = B - N J = D - M
如果组件仍然溢出，可以recursion地将它们分解为更小的组件，然后重新组合。
为了完整起见，由于没有人提到它，所以一些编译器（如GCC）现在实际上为您提供了一个128位的整数。
因此一个简单的解决scheme可能是
(long long)((__int128)A * B - (__int128)C * D)
AB-CD = (AB-CD) * AC / AC = (B/CD/A)*A*C B/C和D/A都不能溢出，所以先计算(B/CD/A) 。 由于最终结果不会根据您的定义溢出，您可以安全地执行剩余的乘法运算并计算(B/CD/A)*A*C ，这是所需的结果。
请注意，如果您的input也可能非常小 ，则B/C或D/A可能会溢出。 如果可能的话，根据input检查可能需要更复杂的操作。
selectK = a big number （例如K = A - sqrt(A) ）
A*B - C*D = (AK)*(BK) - (CK)*(DK) + K*(A-C+BD); // Avoid overflow.
为什么？
(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)
请注意，因为A，B，C和D是大数字，所以AC和BD是小数字。
我可能没有覆盖所有的边缘情况，我也没有严格testing这个，但是这实现了我在80年代尝试在16位cpu上执行32位整数运算时使用的技术。 本质上，你将32位分成两个16位单元并分别使用它们。
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()); } }
打印：
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
这看起来像它的工作。
我敢打赌，我已经错过了一些微妙之处，比如注意符号溢出等，但是我认为它的本质就在那里。
