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是任意自然数。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是任意自然数。我想这似乎太琐碎了。
但是 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
A、B、C、D
是负数,那么E
或F
会更大吗? - Supr最简单、最通用的解决方案是使用一种不能溢出的表示方法,可以通过使用长整型库(例如: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
if( abs( (double)A*B - (double)C*D ) > MAX_LLONG )
Overflow
else
return A*B-C*D;
这种方法的问题在于,您受到双倍精度(54位)的尾数精度的限制,因此需要将乘积A*B和C*D限制在63+54位(或者可能略少一些)。
A
、B
、C
和 D
恰好是相对质数,则它们将没有公共因数)。long double
,取对数似乎是个不错的选择。在这种情况下,可以达到可以接受的精度(并且结果可以四舍五入)。 - user529758123
可以表示为: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 * 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
#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;
}