我需要计算一个看起来像: 的表达式
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, whereMAX = LLONG_MAX - n和 n - 一些自然数。
问问题
8664 次
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
于 2012-11-05T17:30:58.220 回答
68
最简单和最通用的解决方案是使用不会溢出的表示，通过使用长整数库（例如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
于 2012-11-05T17:21:24.503 回答
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溢出行为被标准保证环绕的位置。最后转换回有符号整数是实现定义的部分，但现在几乎可以在所有环境中使用。
如果您需要更多迂腐的解决方案，我认为您必须使用“长算术”
于 2012-11-05T17:57:54.560 回答
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)
于 2012-11-05T17:46:37.640 回答
10
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
于 2012-11-06T06:15:13.700 回答
9
您可以考虑计算所有值的最大公因数，然后在进行算术运算之前将它们除以该因数，然后再次相乘。然而，这假设存在这样的因子（例如，如果A、和恰好是相对素数，则它们将没有公因子）。BCD
同样，您可以考虑使用对数刻度，但这会有点吓人，取决于数值精度。
于 2012-11-05T17:24:56.637 回答
9
如果结果适合 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 位（或者可能更少）。
于 2012-11-05T20:43:37.080 回答
7
您可以将每个数字写在一个数组中，每个元素都是一个数字，并将计算作为polynomials进行。取得到的多项式，它是一个数组，并通过将数组的每个元素乘以 10 到数组中位置的幂来计算结果（第一个位置是最大的，最后一个是零）。
数123可以表示为：
123 = 100 * 1 + 10 * 2 + 3
您只需为其创建一个数组[1 2 3]。
您对所有数字 A、B、C 和 D 执行此操作，然后将它们作为多项式相乘。一旦你得到了多项式，你只需从中重建数字。
于 2012-11-05T17:26:08.487 回答
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_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_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
于 2012-11-05T23:42:53.363 回答
3
如果您知道最终结果可以用您的整数类型表示，则可以使用以下代码快速执行此计算。因为 C 标准规定无符号算术是模算术并且不会溢出，所以您可以使用无符号类型来执行计算。
下面的代码假设有一个相同宽度的无符号类型，并且有符号类型使用所有位模式来表示值（没有陷阱表示，有符号类型的最小值是无符号类型模数一半的负数）。如果这在 C 实现中不成立，则可以为此对 ConvertToSigned 例程进行简单调整。
下面使用signed charandunsigned char来演示代码。对于您的实现，更改to的定义和Signedtotypedef signed long long int Signed;的定义。Unsignedtypedef 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;
}
于 2012-11-05T18:41:32.137 回答
2
您可以尝试将等式分解为不会溢出的较小组件。
AB - CD
= [ A(B - N) - C( D - M )] + [AN - CM]
= ( AK - CJ ) + ( AN - CM)
where K = B - N
J = D - M
如果组件仍然溢出，您可以递归地将它们分解成更小的组件，然后重新组合。
于 2012-11-05T17:51:59.363 回答
2
我可能没有涵盖所有边缘情况，也没有对此进行严格测试，但这实现了我记得在 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
在我看来它正在工作。
我敢打赌我错过了一些微妙之处，例如注意标志溢出等，但我认为本质就在那里。
于 2012-11-07T15:29:19.023 回答
2
为了完整起见，由于没有人提到它，现在一些编译器（例如 GCC）实际上为您提供了一个 128 位整数。
因此，一个简单的解决方案可能是：
(long long)((__int128)A * B - (__int128)C * D)
于 2013-06-28T19:14:31.930 回答
1
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/CorD/A可能会溢出。如果可能，根据输入检查可能需要更复杂的操作。
于 2012-11-05T17:28:03.000 回答
1
选择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是小数。
于 2012-11-07T10:52:45.907 回答