Улавливание и вычисление переполнения при умножении двух больших целых чисел
-
06-07-2019 - |
Вопрос
Я ищу эффективное (необязательно стандартное, элегантное и простое в реализации) решение для умножения относительно больших чисел и сохранения результата в одно или несколько целых чисел :
Допустим, у меня есть два 64-битных целых числа, объявленных следующим образом :
uint64_t a = xxx, b = yyy;
Когда я это делаю a * b
, как я могу определить, приводит ли операция к переполнению, и в этом случае сохранить перенос где-нибудь?
Пожалуйста, обратите внимание, что Я не хочу использовать какую-либо библиотеку с большим количеством номеров поскольку у меня есть ограничения на способ хранения чисел.
Решение
1.Обнаружение переполнения:
x = a * b;
if (a != 0 && x / a != b) {
// overflow handling
}
Редактировать:Фиксированное деление на 0
(спасибо, Марк!)
2.Вычисление переноса это довольно запутанно.Один из подходов заключается в разделении обоих операндов на полслова, а затем применении длинное умножение с полуслова:
uint64_t hi(uint64_t x) {
return x >> 32;
}
uint64_t lo(uint64_t x) {
return ((1L << 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
Чтобы вычислить перенос, мы можем использовать подход, чтобы разбить число на две 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 |
+------------------+----------+----------+----------+----------+----------+
Основываясь на этих результатах, мы можем сделать несколько выводов:
- Очевидно, что подход, основанный на разделении, хотя и прост и переносим, является медленным.
- Ни одна техника не является явным победителем во всех случаях.
- В современных компиляторах подход use-a-large-int является лучшим, если вы можете его использовать
- В старых компиляторах лучше всего использовать подход с длительным умножением
- Удивительно, но GCC 4.9.0 имеет регрессии производительности по сравнению с GCC 4.2.1, а GCC 4.2.1 имеет регрессии производительности по сравнению с GCC 3.3
Версия, которая также работает, когда == 0:
x = a * b;
if (a != 0 && x / a != b) {
// overflow handling
}
Если вам нужно не только обнаружить переполнение, но и захватить перенос, лучше разбить свои числа на 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...
}
Проблема не только в частичных продуктах, но и в том, что любая из сумм может переполниться.
Если бы мне пришлось сделать это по-настоящему, я написал бы подпрограмму расширенного умножения на языке локальной сборки. Например, умножив два 64-разрядных целых числа, чтобы получить 128- битовый результат, который хранится в двух 64-битных регистрах. Все разумные аппаратные средства обеспечивают эту функциональность в единой собственной инструкции умножения & # 8212; она не просто доступна из & Nbsp; C.
Это один из тех редких случаев, когда наиболее элегантным и простым в программировании решением является использование ассемблера. Но это, конечно, не переносимо: - (
Я работал над этой проблемой в эти дни, и я должен сказать, что она произвела на меня впечатление, когда я видел, как люди говорили, что лучший способ узнать, произошло ли переполнение, - это разделить результат, вот и все. совершенно неэффективно и ненужно. Смысл этой функции в том, что она должна быть максимально быстрой. Р>
Существует два варианта обнаружения переполнения:
1 & # 186; - Если возможно, создайте переменную результата в два раза больше, чем множители, например:
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 & # 186; - невозможно создать результирующую переменную, в два раза большую, чем переменная множителей: Тогда вам стоит поиграть с условиями, чтобы определить лучший путь. Продолжая с примером:
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);
}
}
Я надеюсь, что этот код поможет вам создать достаточно эффективную программу, и я надеюсь, что код понятен, если нет, я добавлю некоторые комментарии.
С наилучшими пожеланиями.
Возможно, лучший способ решить эту проблему - это иметь функцию, которая умножает два UInt64 и выдает пару UInt64, верхнюю часть и нижнюю часть результата UInt128. Вот решение, включая функцию, которая отображает результат в шестнадцатеричном виде. Я думаю, вы, возможно, предпочитаете решение C ++, но у меня есть работающее Swift-решение, которое показывает, как решить проблему:
func hex128 (_ hi: UInt64, _ lo: UInt64) -> String
{
var s: String = String(format: "%08X", hi >> 32)
+ String(format: "%08X", hi & 0xFFFFFFFF)
+ String(format: "%08X", lo >> 32)
+ String(format: "%08X", lo & 0xFFFFFFFF)
return (s)
}
func mul64to128 (_ multiplier: UInt64, _ multiplicand : UInt64)
-> (result_hi: UInt64, result_lo: UInt64)
{
let x: UInt64 = multiplier
let x_lo: UInt64 = (x & 0xffffffff)
let x_hi: UInt64 = x >> 32
let y: UInt64 = multiplicand
let y_lo: UInt64 = (y & 0xffffffff)
let y_hi: UInt64 = y >> 32
let mul_lo: UInt64 = (x_lo * y_lo)
let mul_hi: UInt64 = (x_hi * y_lo) + (mul_lo >> 32)
let mul_carry: UInt64 = (x_lo * y_hi) + (mul_hi & 0xffffffff)
let result_hi: UInt64 = (x_hi * y_hi) + (mul_hi >> 32) + (mul_carry >> 32)
let result_lo: UInt64 = (mul_carry << 32) + (mul_lo & 0xffffffff)
return (result_hi, result_lo)
}
Вот пример, чтобы убедиться, что функция работает:
var c: UInt64 = 0
var d: UInt64 = 0
(c, d) = mul64to128(0x1234567890123456, 0x9876543210987654)
// 0AD77D742CE3C72E45FD10D81D28D038 is the result of the above example
print(hex128(c, d))
(c, d) = mul64to128(0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF)
// FFFFFFFFFFFFFFFE0000000000000001 is the result of the above example
print(hex128(c, d))
Вот трюк для определения переполнения умножения двух целых чисел без знака.
Мы отмечаем, что если мы умножим двоичное число шириной N бит на двоичное число шириной M бита, произведение не будет содержать более N + M битов.
Например, если нас попросят умножить трехбитное число на двадцать девять разрядных, мы знаем, что это не переполняет тридцать два бита.
#include <stdlib.h>
#include <stdio.h>
int might_be_mul_oflow(unsigned long a, unsigned long b)
{
if (!a || !b)
return 0;
a = a | (a >> 1) | (a >> 2) | (a >> 4) | (a >> 8) | (a >> 16) | (a >> 32);
b = b | (b >> 1) | (b >> 2) | (b >> 4) | (b >> 8) | (b >> 16) | (b >> 32);
for (;;) {
unsigned long na = a << 1;
if (na <= a)
break;
a = na;
}
return (a & b) ? 1 : 0;
}
int main(int argc, char **argv)
{
unsigned long a, b;
char *endptr;
if (argc < 3) {
printf("supply two unsigned long integers in C form\n");
return EXIT_FAILURE;
}
a = strtoul(argv[1], &endptr, 0);
if (*endptr != 0) {
printf("%s is garbage\n", argv[1]);
return EXIT_FAILURE;
}
b = strtoul(argv[2], &endptr, 0);
if (*endptr != 0) {
printf("%s is garbage\n", argv[2]);
return EXIT_FAILURE;
}
if (might_be_mul_oflow(a, b))
printf("might be multiplication overflow\n");
{
unsigned long c = a * b;
printf("%lu * %lu = %lu\n", a, b, c);
if (a != 0 && c / a != b)
printf("confirmed multiplication overflow\n");
}
return 0;
}
Небольшое количество тестов: (в 64-битной системе):
$ ./uflow 0x3 0x3FFFFFFFFFFFFFFF 3 * 4611686018427387903 = 13835058055282163709 $ ./uflow 0x7 0x3FFFFFFFFFFFFFFF might be multiplication overflow 7 * 4611686018427387903 = 13835058055282163705 confirmed multiplication overflow $ ./uflow 0x4 0x3FFFFFFFFFFFFFFF might be multiplication overflow 4 * 4611686018427387903 = 18446744073709551612 $ ./uflow 0x5 0x3FFFFFFFFFFFFFFF might be multiplication overflow 5 * 4611686018427387903 = 4611686018427387899 confirmed multiplication overflow
Шаги в might_be_mul_oflow
почти наверняка медленнее, чем просто деление теста, по крайней мере на основных процессорах, используемых на настольных рабочих станциях, серверах и мобильных устройствах. На чипах без хорошей поддержки деления это может быть полезно.
Мне приходит в голову, что есть еще один способ сделать этот тест на раннее отклонение.
<Ол> Мы начинаем с пары чисел arng
и brng
, которые инициализируются как 0x7FFF...FFFF
и 1
. Р>
Если a <= arng
и b <= brng
, мы можем заключить, что переполнения нет.
В противном случае мы смещаем 0x3FFF...FFFF
вправо и смещаем 3
влево, добавляя один бит к <=>, чтобы они были <=> и <=>.
Если <=> равно нулю, закончите; в противном случае повторите в 2.
Теперь функция выглядит следующим образом:
int might_be_mul_oflow(unsigned long a, unsigned long b)
{
if (!a || !b)
return 0;
{
unsigned long arng = ULONG_MAX >> 1;
unsigned long brng = 1;
while (arng != 0) {
if (a <= arng && b <= brng)
return 0;
arng >>= 1;
brng <<= 1;
brng |= 1;
}
return 1;
}
}
Если вы просто хотите обнаружить переполнение, как насчет преобразования в double, выполнения умножения и, если
|x| < 2^53, преобразовать в int64
|x| < 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;
}