次のような式を計算する必要があります:
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 - 自然数。
これはあまりにも些細なことだと思います。しかし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 の場合、減算操作中に注意が必要になる場合があります。
たとえば、質問にあるケースでは、1回の繰り返しで済みますが、
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/ ) を使用するか、構造体または配列を使用して表現し、一種の長い乗算 (つまり、各数値を 2 つの 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
結果が long long int に収まる場合、式 A*BC*D は 2^64 を法とする算術演算を実行するので問題なく、正しい結果が得られます。問題は、結果が long long int に収まるかどうかを知ることです。これを検出するには、double を使用した次のトリックを使用できます。
if( abs( (double)A*B - (double)C*D ) > MAX_LLONG )
Overflow
else
return A*B-C*D;
このアプローチの問題は、double の仮数の精度 (54 ビット?) によって制限されるため、A*B と C*D の積を 63+54 ビット (またはおそらくそれ以下) に制限する必要があることです。
すべての値の最大公約数を計算し、算術演算を行う前にそれらをその係数で割り、再度乗算することを検討できます。ただし、これはそのような因数が存在することを前提としています (たとえば、A、B、CおよびDが互いに素である場合、共通の因数はありません)。
同様に、対数スケールで作業することを検討することもできますが、これは数値の精度に左右されるため、少し怖いものになるでしょう。
各数値を配列に書き込み、各要素を数字にして、計算を多項式として行うことができます。結果として得られる多項式 (配列) を取得し、配列の各要素に 10 を配列内の位置で乗じて結果を計算します (最初の位置が最大で、最後の位置が 0 です)。
数は次の123ように表すことができます。
123 = 100 * 1 + 10 * 2 + 3
配列を作成するだけ[1 2 3]です。
これをすべての数値 A、B、C、D に対して行い、それらを多項式として乗算します。結果の多項式が得られたら、そこから数値を再構築するだけです。
asigned long long intは成り立たないがA*B、そのうちの 2 つは成り立つ。したがってA*B、異なる指数のツリー項に分解でき、いずれも 1 に適合し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。
まっすぐな方法に従うと、それぞれに追加のキャリー ビット (正確には 1 ビットの整数) を使用して、 のすべてのペアに対して減算を行うことができAB_iますCD_i。したがって、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の上半分をto に転送することから続けますE_20(32 だけシフトして加算し、上半分を消去しE_10ます)。
E_11これで、(キャリー以外の部分から取得した) 正しい符号を に追加することで、キャリー ビットを取り除くことができますE_20。これがオーバーフローを引き起こした場合、結果も収まりません。
E_10E_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てコードを示します。実装のために、 to の定義と to の定義を変更Signedします。typedef 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;
}
方程式をオーバーフローしない小さなコンポーネントに分割してみることができます。
AB - CD
= [ A(B - N) - C( D - M )] + [AN - CM]
= ( AK - CJ ) + ( AN - CM)
where K = B - N
J = D - M
コンポーネントがまだオーバーフローしている場合は、それらを再帰的に小さなコンポーネントに分割してから再結合できます。
すべてのエッジ ケースをカバーしたわけではなく、これを厳密にテストしたわけでもありませんが、これは、16 ビット CPU で 32 ビット整数演算を実行しようとしたときに 80 年代に使用したことを覚えている手法を実装しています。基本的に、32 ビットを 2 つの 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)
選択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は小さな数であることに注意してください。
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がオーバーフローする可能性があることに注意してください。可能であれば、入力検査に従って、より複雑な操作が必要になる可能性があります。