比較的大きな数を乗算し、結果を1つまたは複数の整数に格納するための効率的な（オプションで標準的でエレガントで実装が簡単な）ソリューションを探しています。
次のように宣言された 2 つの 64 ビット整数があるとします。
uint64_t a = xxx, b = yyy;
私が行うときa * b、操作がオーバーフローになったかどうかをどのように検出し、この場合キャリーをどこかに保存できますか?
数値を格納する方法に制約があるため、多数のライブラリを使用したくないことに注意してください。
1. オーバーフローの検出:
x = a * b;
if (a != 0 && x / a != b) {
// overflow handling
}
編集：除算を修正しました0（マークに感謝します！）
2. キャリーの計算は非常に複雑です。1 つの方法は、両方のオペランドをハーフワードに分割し、長い乗算をハーフワードに適用することです。
uint64_t hi(uint64_t x) {
return x >> 32;
}
uint64_t lo(uint64_t x) {
return ((1ULL << 32) - 1) & x;
}
void multiply(uint64_t a, uint64_t b) {
// actually uint32_t would do, but the casting is annoying
uint64_t s0, s1, s2, s3;
uint64_t x = lo(a) * lo(b);
s0 = lo(x);
x = hi(a) * lo(b) + hi(x);
s1 = lo(x);
s2 = hi(x);
x = s1 + lo(a) * hi(b);
s1 = lo(x);
x = s2 + hi(a) * hi(b) + hi(x);
s2 = lo(x);
s3 = hi(x);
uint64_t result = s1 << 32 | s0;
uint64_t carry = s3 << 32 | s2;
}
部分和自体がオーバーフローしないことを確認するために、最悪のケースを考えます。
x = s2 + hi(a) * hi(b) + hi(x)
しましょうB = 1 << 32。次に、
x <= (B - 1) + (B - 1)(B - 1) + (B - 1)
<= B*B - 1
< B*B
私はこれがうまくいくと信じています - 少なくともSjlverのテストケースを処理します. それを除けば、それはテストされていません (手元に C++ コンパイラがないので、コンパイルすらできないかもしれません)。
アイデアは、積分操作に当てはまる次の事実を使用することです。
a*b > c場合に限りa > c/b
/は整数除算です。
正数のオーバーフローをチェックする疑似コードは次のとおりです。
if (a > max_int64 / b) then "overflow" でなければ "ok" .
ゼロと負の数を処理するには、さらにチェックを追加する必要があります。
非負の C コードはa次のbとおりです。
if (b > 0 && a > 18446744073709551615 / b) {
// overflow handling
}; else {
c = a * b;
}
ノート：
18446744073709551615 == (1<<64)-1
キャリーを計算するには、紙の上でこれを行うように、数値を 2 つの 32 桁に分割し、それらを乗算するアプローチを使用できます。オーバーフローを避けるために数値を分割する必要があります。
コードは次のとおりです。
// split input numbers into 32-bit digits
uint64_t a0 = a & ((1LL<<32)-1);
uint64_t a1 = a >> 32;
uint64_t b0 = b & ((1LL<<32)-1);
uint64_t b1 = b >> 32;
// The following 3 lines of code is to calculate the carry of d1
// (d1 - 32-bit second digit of result, and it can be calculated as d1=d11+d12),
// but to avoid overflow.
// Actually rewriting the following 2 lines:
// uint64_t d1 = (a0 * b0 >> 32) + a1 * b0 + a0 * b1;
// uint64_t c1 = d1 >> 32;
uint64_t d11 = a1 * b0 + (a0 * b0 >> 32);
uint64_t d12 = a0 * b1;
uint64_t c1 = (d11 > 18446744073709551615 - d12) ? 1 : 0;
uint64_t d2 = a1 * b1 + c1;
uint64_t carry = d2; // needed carry stored here
この質問には他にもいくつかの回答がありましたが、そのうちのいくつかには完全にテストされていないコードがあり、これまでのところ、さまざまな可能なオプションを適切に比較した人はいません。
そのため、いくつかの可能な実装を作成してテストしました (最後の実装は、OpenBSD のこのコードに基づいており、Reddit のこちらで説明されています)。コードは次のとおりです。
/* Multiply with overflow checking, emulating clang's builtin function
*
* __builtin_umull_overflow
*
* This code benchmarks five possible schemes for doing so.
*/
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <limits.h>
#ifndef BOOL
#define BOOL int
#endif
// Option 1, check for overflow a wider type
// - Often fastest and the least code, especially on modern compilers
// - When long is a 64-bit int, requires compiler support for 128-bits
// ints (requires GCC >= 3.0 or Clang)
#if LONG_BIT > 32
typedef __uint128_t long_overflow_t ;
#else
typedef uint64_t long_overflow_t;
#endif
BOOL
umull_overflow1(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
long_overflow_t prod = (long_overflow_t)lhs * (long_overflow_t)rhs;
*result = (unsigned long) prod;
return (prod >> LONG_BIT) != 0;
}
// Option 2, perform long multiplication using a smaller type
// - Sometimes the fastest (e.g., when mulitply on longs is a library
// call).
// - Performs at most three multiplies, and sometimes only performs one.
// - Highly portable code; works no matter how many bits unsigned long is
BOOL
umull_overflow2(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
const unsigned long HALFSIZE_MAX = (1ul << LONG_BIT/2) - 1ul;
unsigned long lhs_high = lhs >> LONG_BIT/2;
unsigned long lhs_low = lhs & HALFSIZE_MAX;
unsigned long rhs_high = rhs >> LONG_BIT/2;
unsigned long rhs_low = rhs & HALFSIZE_MAX;
unsigned long bot_bits = lhs_low * rhs_low;
if (!(lhs_high || rhs_high)) {
*result = bot_bits;
return 0;
}
BOOL overflowed = lhs_high && rhs_high;
unsigned long mid_bits1 = lhs_low * rhs_high;
unsigned long mid_bits2 = lhs_high * rhs_low;
*result = bot_bits + ((mid_bits1+mid_bits2) << LONG_BIT/2);
return overflowed || *result < bot_bits
|| (mid_bits1 >> LONG_BIT/2) != 0
|| (mid_bits2 >> LONG_BIT/2) != 0;
}
// Option 3, perform long multiplication using a smaller type (this code is
// very similar to option 2, but calculates overflow using a different but
// equivalent method).
// - Sometimes the fastest (e.g., when mulitply on longs is a library
// call; clang likes this code).
// - Performs at most three multiplies, and sometimes only performs one.
// - Highly portable code; works no matter how many bits unsigned long is
BOOL
umull_overflow3(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
const unsigned long HALFSIZE_MAX = (1ul << LONG_BIT/2) - 1ul;
unsigned long lhs_high = lhs >> LONG_BIT/2;
unsigned long lhs_low = lhs & HALFSIZE_MAX;
unsigned long rhs_high = rhs >> LONG_BIT/2;
unsigned long rhs_low = rhs & HALFSIZE_MAX;
unsigned long lowbits = lhs_low * rhs_low;
if (!(lhs_high || rhs_high)) {
*result = lowbits;
return 0;
}
BOOL overflowed = lhs_high && rhs_high;
unsigned long midbits1 = lhs_low * rhs_high;
unsigned long midbits2 = lhs_high * rhs_low;
unsigned long midbits = midbits1 + midbits2;
overflowed = overflowed || midbits < midbits1 || midbits > HALFSIZE_MAX;
unsigned long product = lowbits + (midbits << LONG_BIT/2);
overflowed = overflowed || product < lowbits;
*result = product;
return overflowed;
}
// Option 4, checks for overflow using division
// - Checks for overflow using division
// - Division is slow, especially if it is a library call
BOOL
umull_overflow4(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
*result = lhs * rhs;
return rhs > 0 && (SIZE_MAX / rhs) < lhs;
}
// Option 5, checks for overflow using division
// - Checks for overflow using division
// - Avoids division when the numbers are "small enough" to trivially
// rule out overflow
// - Division is slow, especially if it is a library call
BOOL
umull_overflow5(unsigned long lhs, unsigned long rhs, unsigned long* result)
{
const unsigned long MUL_NO_OVERFLOW = (1ul << LONG_BIT/2) - 1ul;
*result = lhs * rhs;
return (lhs >= MUL_NO_OVERFLOW || rhs >= MUL_NO_OVERFLOW) &&
rhs > 0 && SIZE_MAX / rhs < lhs;
}
#ifndef umull_overflow
#define umull_overflow2
#endif
/*
* This benchmark code performs a multiply at all bit sizes,
* essentially assuming that sizes are logarithmically distributed.
*/
int main()
{
unsigned long i, j, k;
int count = 0;
unsigned long mult;
unsigned long total = 0;
for (k = 0; k < 0x40000000 / LONG_BIT / LONG_BIT; ++k)
for (i = 0; i != LONG_MAX; i = i*2+1)
for (j = 0; j != LONG_MAX; j = j*2+1) {
count += umull_overflow(i+k, j+k, &mult);
total += mult;
}
printf("%d overflows (total %lu)\n", count, total);
}
これは、私が持っているさまざまなコンパイラとシステムでテストした結果です (この場合、すべてのテストは OS X で行われましたが、結果は BSD または Linux システムでも同様になるはずです)。
+------------------+----------+----------+----------+----------+----------+
| | Option 1 | Option 2 | Option 3 | Option 4 | Option 5 |
| | BigInt | LngMult1 | LngMult2 | Div | OptDiv |
+------------------+----------+----------+----------+----------+----------+
| Clang 3.5 i386 | 1.610 | 3.217 | 3.129 | 4.405 | 4.398 |
| GCC 4.9.0 i386 | 1.488 | 3.469 | 5.853 | 4.704 | 4.712 |
| GCC 4.2.1 i386 | 2.842 | 4.022 | 3.629 | 4.160 | 4.696 |
| GCC 4.2.1 PPC32 | 8.227 | 7.756 | 7.242 | 20.632 | 20.481 |
| GCC 3.3 PPC32 | 5.684 | 9.804 | 11.525 | 21.734 | 22.517 |
+------------------+----------+----------+----------+----------+----------+
| Clang 3.5 x86_64 | 1.584 | 2.472 | 2.449 | 9.246 | 7.280 |
| GCC 4.9 x86_64 | 1.414 | 2.623 | 4.327 | 9.047 | 7.538 |
| GCC 4.2.1 x86_64 | 2.143 | 2.618 | 2.750 | 9.510 | 7.389 |
| GCC 4.2.1 PPC64 | 13.178 | 8.994 | 8.567 | 37.504 | 29.851 |
+------------------+----------+----------+----------+----------+----------+
これらの結果に基づいて、いくつかの結論を導き出すことができます。
== 0 の場合にも機能するバージョン:
x = a * b;
if (a != 0 && x / a != b) {
// overflow handling
}
clang と gcc で簡単かつ高速:
unsigned long long t a, b, result;
if (__builtin_umulll_overflow(a, b, &result)) {
// overflow!!
}
これにより、利用可能な場合はオーバーフロー検出のハードウェア サポートが使用されます。コンパイラの拡張機能であるため、符号付き整数のオーバーフロー (umul を smul に置き換える) を処理することもできますが、これは C++ では未定義の動作です。
オーバーフローを検出するだけでなく、キャリーをキャプチャする必要がある場合は、数値を 32 ビットの部分に分割することをお勧めします。コードは悪夢です。以下は単なるスケッチです。
#include <stdint.h>
uint64_t mul(uint64_t a, uint64_t b) {
uint32_t ah = a >> 32;
uint32_t al = a; // truncates: now a = al + 2**32 * ah
uint32_t bh = b >> 32;
uint32_t bl = b; // truncates: now b = bl + 2**32 * bh
// a * b = 2**64 * ah * bh + 2**32 * (ah * bl + bh * al) + al * bl
uint64_t partial = (uint64_t) al * (uint64_t) bl;
uint64_t mid1 = (uint64_t) ah * (uint64_t) bl;
uint64_t mid2 = (uint64_t) al * (uint64_t) bh;
uint64_t carry = (uint64_t) ah * (uint64_t) bh;
// add high parts of mid1 and mid2 to carry
// add low parts of mid1 and mid2 to partial, carrying
// any carry bits into carry...
}
問題は部分積だけではなく、いずれかの合計がオーバーフローする可能性があるという事実です。
これを実際に行う必要がある場合は、拡張乗算ルーチンをローカル アセンブリ言語で記述します。 つまり、たとえば、2 つの 64 ビット整数を乗算して 128 ビットの結果を取得し、2 つの 64 ビット レジスタに格納します。すべての適切なハードウェアは、単一のネイティブ乗算命令でこの機能を提供します。C からだけアクセスできるわけではありません。
これは、実際にアセンブリ言語を使用することが最もエレガントでプログラミングが簡単な解決策である、まれなケースの 1 つです。しかし、それは確かに移植可能ではありません:-(
私は今日この問題に取り組んできましたが、オーバーフローがあったかどうかを知る最善の方法は結果を分割することであると人々が言っているのを何度も見たことに感銘を受けました。それは完全に非効率的であり、不要。この関数のポイントは、可能な限り高速でなければならないということです。
オーバーフロー検出には 2 つのオプションがあります。
1º- 可能であれば、乗数の 2 倍の大きさの結果変数を作成します。次に例を示します。
struct INT32struct {INT16 high, low;};
typedef union
{
struct INT32struct s;
INT32 ll;
} INT32union;
INT16 mulFunction(INT16 a, INT16 b)
{
INT32union result.ll = a * b; //32Bits result
if(result.s.high > 0)
Overflow();
return (result.s.low)
}
オーバーフローが発生したかどうかはすぐにわかり、コードは機械語で書かなくても可能な限り高速です。コンパイラによっては、このコードをマシン コードで改善できます。
2º- 乗数変数の 2 倍の大きさの結果変数を作成することは不可能です。最適なパスを決定するには、if 条件を使用する必要があります。例を続ける：
INT32 mulFunction(INT32 a, INT32 b)
{
INT32union s_a.ll = abs(a);
INT32union s_b.ll = abs(b); //32Bits result
INT32union result;
if(s_a.s.hi > 0 && s_b.s.hi > 0)
{
Overflow();
}
else if (s_a.s.hi > 0)
{
INT32union res1.ll = s_a.s.hi * s_b.s.lo;
INT32union res2.ll = s_a.s.lo * s_b.s.lo;
if (res1.hi == 0)
{
result.s.lo = res1.s.lo + res2.s.hi;
if (result.s.hi == 0)
{
result.s.ll = result.s.lo << 16 + res2.s.lo;
if ((a.s.hi >> 15) ^ (b.s.hi >> 15) == 1)
{
result.s.ll = -result.s.ll;
}
return result.s.ll
}else
{
Overflow();
}
}else
{
Overflow();
}
}else if (s_b.s.hi > 0)
{
//Same code changing a with b
}else
{
return (s_a.lo * s_b.lo);
}
}
このコードが非常に効率的なプログラムを作成するのに役立つことを願っています。また、コードが明確であることを願っています。そうでない場合は、コメントを追加します。
よろしくお願いします。
オーバーフローを検出したいだけなら、double に変換して乗算を行い、
|×| < 2^53、int64 に変換
|×| < 2^63、int64 を使用して乗算を行います
それ以外の場合は、必要なエラーが発生しますか?
これはうまくいくようです：
int64_t safemult(int64_t a, int64_t b) {
double dx;
dx = (double)a * (double)b;
if ( fabs(dx) < (double)9007199254740992 )
return (int64_t)dx;
if ( (double)INT64_MAX < fabs(dx) )
return INT64_MAX;
return a*b;
}