如何在 expr 中避免溢出

最简单也是最通用的解决方案是使用一个不会溢出的表示，或者使用一个长整数库(例如 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位

您可以考虑为所有值计算一个最大的公因数，然后在进行算术运算之前将它们除以该公因数，然后再次进行乘法运算。但是，这里假设存在这样一个因子(例如，如果 A、 B、 C和 D碰巧是相对质数，它们就不会有一个公因子)。

类似地，您可以考虑使用对数刻度，但这将有点可怕，这取决于数值的精度。

您可以将每个数字写入一个数组中，每个元素为一个数字，并以 多项式的形式进行计算。取得的多项式是一个数组，通过将数组中的每个元素乘以10乘以数组中位置的幂(第一个位置是最大的，最后一个位置是零)来计算结果。

数字 123可以表示为:

123 = 100 * 1 + 10 * 2 + 3

为其创建一个数组 [1 2 3]。

对所有数字 A、 B、 C 和 D 都这样做，然后将它们作为多项式相乘。一旦你得到了结果多项式，你只需要从中重建数字。

AB-CD = (AB-CD) * AC / AC = (B/C-D/A)*A*C.B/C和 D/A都不会溢出，所以先计算 (B/C-D/A)。由于最终结果不会根据您的定义溢出，因此您可以安全地执行剩余的乘法并计算 (B/C-D/A)*A*C，这是所需的结果。

注意，如果您的输入也可以是 非常小，则 B/C或 D/A可能会溢出。如果可能的话，根据输入检查，可能需要更复杂的操作。

我想这似乎太琐碎了。 但 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 所指出的，如果类型是 unsignedlong long，那么在减法操作期间可能需要注意。

例如，对于问题中的情况，只需要一次迭代,

 MAX * MAX - (MAX - 1) * (MAX + 1)
A     B       C           D


E = B - D = -1
F = C - A = -1


AE - DF = {MAX * -1} - {(MAX + 1) * -1} = -MAX + MAX + 1 = 1

这应该会奏效(我认为) :

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)

你可以试着把这个方程分解成不会溢出的小分量。

AB - CD
= [ A(B - N) - C( D - M )] + [AN - CM]


= ( AK - CJ ) + ( AN - CM)


where K = B - N
J = D - M

如果组件仍然溢出，您可以递归地将它们分解成更小的组件，然后重新组合。

请注意，这并不是标准的，因为它依赖于环绕式签名溢出(GCC 有启用此功能的编译器标志)

但是如果你只是在 long long中做所有的计算，直接应用公式的结果是: < br > (A * B - C * D)将是准确的，只要正确的结果符合 long long。

这里有一个解决方案，它只依赖于实现定义的将无符号整数转换为有符号整数的行为。但是现在几乎所有的系统都可以使用这种方法。

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

这将输入强制转换为 unsigned long long，其中溢出行为保证由标准包装。最后转换回一个有符号整数是实现定义的部分，但是今天几乎在所有环境中都可以工作。

如果你需要更迂腐的解决方案，我认为你必须使用“长算术”

如果 知道最终结果可以用整数类型表示，则可以使用下面的代码快速执行此计算。因为 C 标准规定无符号算术是模算术，不会溢出，所以您可以使用无符号类型来执行计算。

下面的代码假设有一个相同宽度的无符号类型，并且有符号类型使用所有位模式来表示值(没有陷阱表示，有符号类型的最小值是无符号类型模的一半的负数)。如果这在 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;
}

如果结果符合一个长长的 int，那么表达式 A * B-C * D 就没问题，因为它执行算术 mod 2 ^ 64，并将给出正确的结果。问题是要知道结果是否适合长整型。为了检测这一点，您可以使用以下技巧使用双精度:

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

这种方法的问题在于，您受到双精度尾数(54位?)的限制因此，您需要将产品 A * B 和 C * D 限制为63 + 54位(或者可能少一点)。

虽然 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同样进行减法，对每对 AB_i和 CD_i使用额外的进位(准确地说是1位整数)。所以如果我们说 E = A * B-C * 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_20(移动32并添加，然后擦除 E_10的上半部分)。

现在您可以通过将带有正确符号(从非进位部分获得)的进位 E_11添加到 E_20来去除它。如果这触发了溢出，结果也不符合。

E_10现在有足够的“空间”从 E_00(移位、加法、擦除)和进位 E_01取出上半部分。

E_10现在可能再次变大，所以我们重复转移到 E_20。

此时，E_20必须变为零，否则结果将不匹配。作为转移的结果，E_10的上半部分也是空的。

最后一步是再次将 E_20的下半部分转移到 E_10。

如果 E=A*B+C*D符合 signed long long int的预期成立，我们现在有

E_20=0
E_10=0
E_00=E

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

选择 K = a big number(例如 K = A - sqrt(A))

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

为什么？

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

注意，因为 A、 B、 C 和 D 是大数，所以 A-C和 B-D是小数。

我可能没有涵盖所有的边缘情况，也没有严格测试过这个，但是这个实现了一种技术，我记得在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 ) {
/*
* e.g. 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 * b.l;
long hbh = h * b.h;
long lbh = l * b.h;
long hbl = h * b.l;
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

在我看来起作用了。

我打赌我已经错过了一些微妙的东西，比如观察符号溢出等，但我认为本质是存在的。

为了完整起见，由于没有人提到它，现在一些编译器(例如 GCC)实际上提供了一个128位整数。

因此，一个简单的解决办法可以是:

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