如何避免 expr 中的溢出。 A B C D

我猜这似乎太微不足道了。但 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

---

最简单和最通用的解决方案是使用不会溢出的表示，通过使用长整数库（例如 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)

这会将输入转换为 unsigned long long，其中溢出行为被标准保证为环绕。最后转换回有符号整数是实现定义的部分，但现在几乎可以在所有环境中使用。

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

---

这应该有效（我认为）：

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)

---

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

---

您可以考虑计算所有值的最大公因数，然后在进行算术运算之前将它们除以该因数，然后再次相乘。然而，这假设存在这样的因数（例如，如果 A、B、C 和 D 恰好是互质数，则它们将没有公因数）。

同样，您可以考虑使用对数刻度，但这会有点吓人，取决于数值精度。

---

如果结果适合 long long int，则表达式 A*BC*D 是可以的，因为它执行算术 mod 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 位（或者可能更少）。

---

您可以将每个数字写入一个数组，每个元素都是一个数字，然后按照 polynomials 进行计算。取得到的多项式，它是一个数组，并通过将数组的每个元素乘以 10 到数组中位置的幂来计算结果（第一个位置是最大的，最后一个是零）。

数 123 可以表示为：

123 = 100 * 1 + 10 * 2 + 3

您只需为其创建一个数组 [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_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

---

如果您知道最终结果可以用您的整数类型表示，则可以使用以下代码快速执行此计算。因为 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;
}

---

您可以尝试将等式分解为不会溢出的较小组件。

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

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

    where K = B - N
          J = D - M

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

---

我可能没有涵盖所有边缘情况，也没有对此进行严格测试，但这实现了我记得在 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)

---

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 可能会溢出。如果可能，根据输入检查可能需要更复杂的操作。

---

选择 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 是小数。

---