如何避免表达式A * B - C * D中的溢出。

161
我需要计算一个表达式,其形式如下: A*B - C*D,其中它们的类型为:signed long long int A, B, C, D; 每个数字都可能非常大(不会溢出其类型)。虽然A*B可能会导致溢出,但同时表达式A*B - C*D可以非常小。那么我该如何正确地计算它呢?
例如:MAX * MAX - (MAX - 1) * (MAX + 1) == 1,其中MAX = LLONG_MAX - n,n是任意自然数。

17
准确性有多重要? - Anirudh Ramanathan
1
@Cthulhu,非常棒的问题。他可以尝试使用更小的数字通过将它们全部除以10或其他数字来制作一个等效的函数,然后再乘以结果。 - Chris
4
变量A、B、C、D是有符号的。这意味着"A - C"可能会发生溢出。这是否需要考虑,或者你是否确定在你的数据中不会发生这种情况? - William Morris
2
@MooingDuck 但是您可以事先检查操作是否会溢出 https://dev59.com/bE7Sa4cB1Zd3GeqP8f-O#3224630 - bradgonesurfing
1
@Chris:不,我的意思是没有一种可移植的方法来检查是否发生了有符号溢出。(Brad 正确指出您可以可移植地检测到它将会发生)。使用内联汇编是众多非可移植检查方式之一。 - Mooing Duck
显示剩余7条评论
15个回答

120

我想这似乎太琐碎了。 但是 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所指出的,如果类型是unsigned long 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

4
@Caleb,将相同的算法应用于“C*D”。 - Chris
2
我认为你应该解释一下E代表什么。 - Caleb
7
长长整型和双精度浮点型都是64位。由于双精度浮点型需要分配一些位数给指数部分,所以它具有比长长整型更小的可能值范围,但不会失去精度。 - Jim Garrison
3
似乎只有在所有数字都非常大的情况下,这种方法才能奏效...例如,当{A,B,C,D} = {MAX,MAX,MAX,2}时,仍然会发生溢出。问题陈述中指出“每个数字可以非常大”,但并不清楚是否必须如此。 - Kevin K
4
如果任何一个值A、B、C、D是负数,那么EF会更大吗? - Supr
显示剩余10条评论

69

最简单、最通用的解决方案是使用一种不能溢出的表示方法,可以通过使用长整型库(例如: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的任何比特。


9
上面的计算并没有看起来那么复杂,其实它只是在基于2^32的简单长乘法运算,使用C语言编写代码应该更加简单。此外,在程序中创建通用函数来完成这项工作是一个很好的主意。 - Ofir

46

请注意,这不是标准方式,因为它依赖于环绕有符号溢出。 (GCC有启用此功能的编译器标志。)

但如果您只是在long long中进行所有计算,则直接应用公式的结果:
(A * B - C * D)只要正确结果适合long long,就会准确。


这里有一个解决方法,仅依赖于将无符号整数强制转换为有符号整数的实现定义行为。 但可以预期这在今天几乎所有系统上都能正常工作。

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

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


如果您需要更严谨的解决方案,我认为您必须使用 "长算术运算"。


+1 你是唯一一个注意到这个问题的人。唯一棘手的部分是设置编译器进行环绕溢出检查,以及检查正确的结果是否适合于“long long”。 - Mysticial
2
即使是最朴素的版本,没有任何技巧,在大多数实现上也能正确执行;虽然这不是标准保证的,但你必须找到一个1's补码机器或其他相当奇怪的设备才能使它失败。 - hobbs
1
我认为这是一个重要的答案。我同意假设实现特定行为可能不是正确的编程方式,但每个工程师都应该了解模数算术以及如何获得正确的编译器标志以确保在性能至关重要时保持一致的行为。DSP工程师依赖于固定点滤波器实现的这种行为,而接受的答案将具有不可接受的性能。 - Peter M

18

这应该可以运作(我想):

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)

1
感谢@bradgonesurfing - 你能提供这样的输入吗?我已经更新了我的答案,执行了它,它可以工作(bd和ca都是0)... - paquetp
1
嗯,现在我想想可能不用。当d = 1且a = 1,b = maxint和c = maxint时,它仍然有效。很酷 :) - bradgonesurfing
3
@MooingDuck,但是你集合的最终答案也溢出了,因此它不是一个有效的设置。只有在每一边都是相同符号时才能工作,这样得到的减法结果才在范围内。 - bradgonesurfing
1
当这个回答相比得分最高的回答分数如此之低时,StackOverflow上有一些奇怪的地方。 - bradgonesurfing
1
这个答案使用了64位除法,效率不是很高。相比之下,Ofir只使用加法和乘法。我预计在各种硬件上,从低端ARM9到顶级x86服务器CPU上,这将比64位除法高一个数量级的效率。 - MSalters
显示剩余3条评论

11
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

9
如果结果适合长长整型(long long int),则表达式A*B-C*D是正确的,因为它执行2^64取模的算术运算,并将给出正确的结果。问题在于如何知道结果是否适合长长整型(long long int)。为了检测这一点,您可以使用以下使用双精度浮点数的技巧:
if( abs( (double)A*B - (double)C*D ) > MAX_LLONG ) 
    Overflow
else 
    return A*B-C*D;

这种方法的问题在于,您受到双倍精度(54位)的尾数精度的限制,因此需要将乘积A*B和C*D限制在63+54位(或者可能略少一些)。


这是最实用的例子。清晰明了,能够给出正确的答案(或在输入有误时抛出异常)。 - Mark Lakata
1
好漂亮!你没有像其他人一样掉进陷阱。只有一件事:我敢打赌,有些例子中双倍计算之下的值仅仅是由于舍入误差而低于MAX_LLONG。我的数学直觉告诉我,你应该计算双倍和长整型结果之间的差异,然后将其与MAX_LLONG/2或其他东西进行比较。这个差异是双倍计算的舍入误差加上溢出,通常应该相对较低,但在我提到的情况下会很大。但现在我太懒了,不想确定。 :-) - Hans-Peter Störr

9
您可以考虑计算所有值的最大公因数,然后在进行算术运算之前将它们除以该因数,然后再进行乘法运算。这假设存在这样一个因数(例如,如果 ABCD 恰好是相对质数,则它们将没有公共因数)。
同样,您可以考虑使用对数比例尺,但这可能会有一些数字精度问题,需要谨慎使用。

1
如果有long double,取对数似乎是个不错的选择。在这种情况下,可以达到可以接受的精度(并且结果可以四舍五入)。 - user529758

7
你可以将每个数字写成一个数组,每个元素都是一个数字,并按照多项式的方式进行计算。取得结果的多项式,即为一个数组,通过将数组中每个元素乘以10的幂次方(第一个位置为最大值,最后一个位置为0)来计算结果。
数字123可以表示为:
123 = 100 * 1 + 10 * 2 + 3

对于这个问题,你只需要创建一个数组[1 2 3]

对于所有的数字A、B、C和D,你都要这样做,然后将它们作为多项式相乘。一旦你得到了结果多项式,你就可以从中重构出数字。


2
不知道那是什么,但我会去找一下。这只是我在和女友逛街时脑海中的一个解决方案 :)。 - Mihai
你正在使用基于十进制数组实现大数。GMP是一个优质的大数库,它使用基于4294967296的底数。速度要快得多。虽然答案正确且有用,但不会点踩。 - Mooing Duck
谢谢 :)。知道这种方法是有用的,但还有更好的方法,所以不要像这样做。至少在这种情况下不要这样做 :)。 - Mihai
无论如何...使用这个解决方案,您可以计算比任何原始类型都要大得多的数字(如100位数),并将结果保留为数组。这值得点赞 :p - Mihai
我不确定它是否会得到赞同,因为这种方法(虽然有效且相对易于理解)需要大量内存并且速度较慢。 - Mooing Duck
让我们在聊天中继续这个讨论:http://chat.stackoverflow.com/rooms/19108/discussion-between-mihai-and-mooing-duck - Mihai

6

虽然一个 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_iCD_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_20中来摆脱进位位E_11。如果这触发了溢出,则结果也不适合。 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

1
这实际上是简化公式,如果使用Ofir的乘法公式并删除每个无用的临时结果,就会得到这个公式。 - dronus

3
如果你知道最终结果可以在整数类型中表示,你可以使用以下代码快速进行计算。因为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;
}

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接