Поймать и вычислить переполнение при умножении двух больших целых чисел
Я ищу эффективное (необязательно стандартное, элегантное и простое в реализации) решение для умножения относительно больших чисел и сохранения результата в одно или несколько целых чисел:
Допустим, у меня есть два 64-битных целых числа, объявленных так:
uint64_t a = xxx, b = yyy;
Когда я делаю a * bКак я могу обнаружить, если операция приводит к переполнению и в этом случае сохранить перенос где-нибудь?
Обратите внимание, что я не хочу использовать какую-либо библиотеку больших чисел, поскольку у меня есть ограничения на способ хранения номеров.
15 ответов
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
Я верю, что это сработает - по крайней мере, он справится с контрольным тестом Сливера. Кроме того, он не проверен (и может даже не скомпилироваться, так как у меня больше нет компилятора C++).
Идея состоит в том, чтобы использовать следующий факт, который справедлив для интегральной операции:
a*b > c если и только если a > c/b
/ является неотъемлемым делением здесь.
Ниже приводится псевдокод для проверки на переполнение для положительных чисел:
if (a> max_int64 / b) тогда "переполнение", иначе "ок".
Для обработки нулей и отрицательных чисел вы должны добавить больше проверок.
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 |
+------------------+----------+----------+----------+----------+----------+
На основании этих результатов можно сделать несколько выводов:
- Очевидно, что подход, основанный на разделении, хотя и прост и переносим, медленен.
- Ни одна техника не является явным победителем во всех случаях.
- На современных компиляторах лучше использовать подход "больше-больше-int", если вы можете его использовать
- На старых компиляторах лучше всего подходит метод длинного умножения
- Удивительно, но GCC 4.9.0 имеет регрессии производительности по GCC 4.2.1, а GCC 4.2.1 имеет регрессии производительности по GCC 3.3
Легко и быстро с clang и gcc:
unsigned long long t a, b, result;
if (__builtin_umulll_overflow(a, b, &result)) {
// overflow!!
}
При этом будет использоваться аппаратная поддержка для обнаружения переполнения, если это возможно. Будучи расширением компилятора, он может даже обрабатывать подписанное целочисленное переполнение (замените umul на smul), хотя это поведение undefined в C++.
Версия, которая также работает, когда == 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-разрядных регистрах. Все разумные аппаратные средства обеспечивают эту функциональность в одной собственной инструкции умножения - она не просто доступна из C.
Это один из тех редких случаев, когда наиболее элегантным и простым в программировании решением является использование ассемблера. Но это конечно не портативно:-(
Библиотека переносимости GNU (Gnulib) содержит модуль intprops, в котором есть макросы, которые эффективно проверяют, будут ли арифметические операции переполнены.
Например, если произойдет переполнение при умножении, INT_MULTIPLY_OVERFLOW (a, b) даст 1.
Возможно, лучший способ решить эту проблему - это иметь функцию, которая умножает два 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 бит или выше.
В общем, вам нужно проверить, достаточно ли велика сумма начальных нулевых битов в обоих факторах, чтобы предотвратить переполнение. Что мне действительно нравится в решении, так это математический аспект умножения битов. Пусть оба операнда и результат равны n бит, тогда легко доказать, что любая сумма ведущих нулей, меньшая, чем разрядная ширина результата n, может (отредактировано благодаря комментарию) дать вам переполнение или, скорее, что оно не переполнится, если ведущий ноль сумма достаточно велика (не менее n). Если вы возьмете наименьшее возможное произведение факторов, сумма ведущих нулей которых вместе сокращает n -1, то переполнение обязательно. (Однако обратите внимание: если начальная нулевая сумма равна n-1, это неоднозначно. Потребуется немного расширить ее, чтобы распознать особый случай начальной нулевой суммы == n-1.)
Причина, по которой этот подход значительно более эффективен, чем предложенный выше метод деления, основана на том факте, что многие популярные процессоры поддерживают подсчет ведущих нулей с помощью собственной машинной инструкции, которая намного быстрее, чем проверка ветвлений для нуля, деление и затем сравнение-ветвление. еще раз. Вычисление делений обычно занимает гораздо больше времени (представьте себе, что на ARM Cortex-M требуется более 10 тактов), чем подсчет ведущих нулей (который выполняется за один цикл, если поддерживается как машинная инструкция), и их подсчет даже не требуется. программируется самостоятельно, даже не с помощью встроенного ассемблера!
Уловка заключается в использовании встроенных / встроенных функций. В GCC это выглядит так:
/**@fn static inline _Bool chk_mul_ov(uint32_t f1, uint32_t f2)
* @return one, if a 32-Bit-overflow occurs when unsigned-unsigned-multipliying f1 with f2 otherwise zero. */
static inline _Bool chk_mul_ov(uint32_t f1, uint32_t f2) {
int lzsum = builtin_clz(f1) + builtin_clz(f2); //leading zero sum
return
lzsum < sizeof(f1)*8-1 || ( //if too small, overflow guaranteed
lzsum == sizeof(f1)*8-1 && //if special case, do further check
(int32_t)((f1 >> 1)*f2 + (f1 & 1)*(f2 >> 1)) < 0 //check product rightshifted by one
);
}
...
if (chk_mul_ov(f1, f2)) {
//error handling
}
...
Это просто пример, который работает только для 32-битного умножения без знака и без знака. И даже не требуется многократного сдвига бит (потому что некоторые микроконтроллеры реализуют только однобитовый сдвиг). Однако, если нет инструкции подсчета-начальных нулей, но существует инструкция умножения, это может быть не лучше, чем просто умножение всех битов.
У других компиляторов есть свой собственный способ определения встроенных функций для операций CLZ. Я думаю, что во многих случаях это решение может быть столь же быстрым, как вычисление 128-битного произведения и проверка его большей половины. Многие процессоры все еще могут даже не предоставлять инструкции типа UUMULL, SUMULL или SSMULL или 64-битные регистры, что означает, что им потребуется 4 регистра для 128-битных. Следовательно, метод подсчета с ведущим нулем может даже лучше масштабироваться (в худшем случае), чем использование высокооптимизированного 128-битного умножения для проверки 64-битного переполнения. Умножение требует линейных накладных расходов, в то время как для подсчета битов нужны только линейные накладные расходы. Я не пробовал свою идею на практике, но надеюсь, что это поможет решить проблему.
Я работал с этой проблемой в эти дни, и я должен сказать, что это произвело на меня впечатление, когда я видел, как люди говорили, что лучший способ узнать, было ли переполнение, разделить результат, что совершенно неэффективно и ненужным. Смысл этой функции в том, что она должна быть максимально быстрой.
Существует два варианта обнаружения переполнения:
1º - Если возможно, создайте переменную результата в два раза больше, чем множители, например:
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º- Невозможно создать результирующую переменную, в два раза большую, чем переменная множителей: тогда вам следует поиграть с условиями для определения наилучшего пути. Продолжая с примером:
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);
}
}
Я надеюсь, что этот код поможет вам получить достаточно эффективную программу, и я надеюсь, что код понятен, если нет, я добавлю некоторые комментарии.
с уважением.
Я думал об этом сегодня и наткнулся на этот вопрос, мои мысли привели меня к такому результату. TL; DR, хотя я нахожу его "элегантным" в том смысле, что он использует только несколько строк кода (легко может быть одним лайнером) и имеет некоторую мягкую математику, которая концептуально упрощает что-то относительно простое, это в основном "интересно", и я не знаю Не тестировал.
Если вы думаете о беззнаковом целом числе как о единственной цифре с основанием 2 ^ n, где n - количество бит в целом числе, вы можете сопоставить эти числа с радианами по единичной окружности, например
radians(x) = x * (2 * pi * rad / 2^n)
Когда целое число переполняется, это эквивалентно обтеканию круга. Таким образом, вычисление переноса эквивалентно вычислению того, сколько раз умножение будет повторяться по кругу. Чтобы вычислить, сколько раз мы оборачиваем круг, мы делим радианы (x) на 2pi радиана. например
wrap(x) = radians(x) / (2*pi*rad)
= (x * (2*pi*rad / 2^n)) / (2*pi*rad / 1)
= (x * (2*pi*rad / 2^n)) * (1 / 2*pi*rad)
= x * 1 / 2^n
= x / 2^n
Что упрощает
wrap(x) = x / 2^n
Это имеет смысл. Сколько раз оборачивается число, например 15 с основанием 10, равно
15 / 10 = 1.5, или в полтора раза. Однако мы не можем использовать здесь 2 цифры (при условии, что мы ограничены одной цифрой 2^64).
Скажем, у нас есть a * b, с основанием R, мы можем вычислить перенос с помощью
Consider that: wrap(a * b) = a * wrap(b)
wrap(a * b) = (a * b) / R
a * wrap(b) = a * (b / R)
a * (b / R) = (a * b) / R
carry = floor(a * wrap(b))
Возьмем к примеру
a = 9 а также
b = 5, которые являются множителями 45 (т. е.
9 * 5 = 45).
wrap(5) = 5 / 10 = 0.5
a * wrap(5) = 9 * 0.5 = 4.5
carry = floor(9 * wrap(5)) = floor(4.5) = 4
Обратите внимание, что если бы перенос был 0, то у нас не было бы переполнения, например, если бы
a = 2,
b=2.
В C / C++ (если компилятор и архитектура это поддерживают) мы должны использовать long double.
Таким образом, мы имеем:
long double wrap = b / 18446744073709551616.0L; // this is b / 2^64
unsigned long carry = (unsigned long)(a * wrap); // floor(a * wrap(b))
bool overflow = carry > 0;
unsigned long c = a * b;
c - это младшая значащая "цифра", т.е. по основанию 10
9 * 9 = 81,
carry = 8, а также
c = 1.
Мне это было интересно в теории, поэтому я подумал, что поделюсь им, но одна серьезная оговорка связана с точностью вычислений с плавающей запятой в компьютерах. Использование long double может привести к ошибкам округления для некоторых чисел при вычислении
wrapзависит от того, сколько значащих цифр использует ваш компилятор / арка для длинных удвоений, я считаю, что для уверенности должно быть на 20 больше. Другая проблема с этим результатом заключается в том, что он может не работать так же хорошо, как некоторые другие решения, просто используя плавающую точку и деление.
Вот трюк для определения переполнения умножения двух целых чисел без знака.
Мы отмечаем, что если мы умножим двоичное число шириной 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 может быть переполнение умножения 7 * 4611686018427387903 = 13835058055282163705 подтвержденное переполнение умножения $ ./uflow 0x4 0x3FFFFFFFFFFFFFFF может быть переполнение умножения 4 * 4611686018427387903 = 18446744073709551612 $./uflow 0x5 0x3FFFFFFFFFFFFFFF может быть переполнение умножения 5 * 4611686018427387903 = 4611686018427387899 подтвержденное переполнение умножения
Шаги в might_be_mul_oflow почти наверняка медленнее, чем просто деление теста, по крайней мере на основных процессорах, используемых в настольных рабочих станциях, серверах и мобильных устройствах. На чипах без хорошей поддержки деления это может быть полезно.
Мне приходит в голову, что есть еще один способ сделать этот тест на раннее отклонение.
Начнем с пары чисел
arngа такжеbrngкоторые инициализированы в0x7FFF...FFFFа также1,Если
a <= arngа такжеb <= brngмы можем сделать вывод, что переполнения нет.В противном случае мы сдвигаемся
arngвправо и сдвигbrngслева, добавив один бит кbrngтак, чтобы они0x3FFF...FFFFа также3,Если
arngноль, финиш; в противном случае повторите в 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;
}
}
Когда вы используете, например, 64-битные переменные, реализуйте «количество значащих битов» с помощью nsb(var) = {64 - clz(var); }.
clz(var) = подсчитывать ведущие нули в var, встроенной команде для GCC и Clang, или, возможно, доступной с встроенной сборкой для вашего процессора.
Теперь используйте тот факт, что nsb(a * b) <= nsb(a) + nsb(b), чтобы проверить переполнение. Чем меньше - всегда на 1 меньше.
Ссылка GCC: встроенная функция: int __builtin_clz (unsigned int x) Возвращает количество ведущих 0-битов в x, начиная с позиции самого старшего бита. Если x равен 0, результат не определен.
на основе ответа Меритона :
использует:
4
2
3
вместо:
8>>
8&
4+
2<<
2|
хотя я не уверен, как проверить "правильность" моего решения
__Bool hasCarryимеет возможность короткого замыкания (если вам не нужен фактический перенос)
uint64_t multiply_compute_carry(uint64_t a, uint64_t b) { //https://stackoverflow.com/questions/1815367/catch-and-compute-overflow-during-multiplication-of-two-large-integers#1815371
uint64_t highA = a >> 32;
uint64_t highB = b >> 32;
uint64_t lowA = a & 0b11111111111111111111111111111111;
uint64_t lowB = b & 0b11111111111111111111111111111111;
// long multiplication of 12*34: 1,2,3,4 are symbols, 001 is padded 1, 002 is padded 2,...
// | 001002
// | 003004
// |
// | 2*4
// | 1*4
// | 2*3
// |1*3
// row1:lowA*lowB
// row2:(highA*lowB)<<32
// row3:(highB*lowA)<<32
// row4:(highA*highB)<<64
// sum:(lowA*lowB) + ((highA*lowB)<<32) + ((highB*lowA)<<32) + ((highA*highB)<<64)
// carry:((lowA*lowB) + ((highA*lowB)<<32) + ((highB*lowA)<<32) + ((highA*highB)<<64)) >> 64
// carry:((lowA*lowB) + ((highA*lowB)<<32) + ((highB*lowA)<<32)) >> 64 + (highA*highB)
// carry:(((lowA*lowB)>>32) + (highA*lowB) + (highB*lowA)) >> 32 + (highA*highB)
// uint64_t carry = highA * highB + ((((lowB*highA) + (lowA*highB) + ((lowA*lowB) >> 32))) >> 32);
// _Bool hasCarry = highA * highB || ((((lowB*highA) + (lowA*highB) + ((lowA*lowB) >> 32))) >> 32);
return highA * highB + ((((lowB*highA) + (lowA*highB) + ((lowA*lowB) >> 32))) >> 32);
}
на основе ответа ChrisoLosoph :
все, что я сделал, это изменил порядок ветвления и предоставил альтернативное объяснение, используя диапазоны
_Bool check_multiplication_overflow(uint64_t a, uint64_t b) { //https://stackoverflow.com/questions/1815367/catch-and-compute-overflow-during-multiplication-of-two-large-integers#59109502
// assume a>0, b>0:
// assume integer log2:
// if log2(a) = A
// then a ∈ [2^A,2^(A+1))
// let log2(a) = A, log2(b) = B
// ab ∈ [2^A*2^B,2^(A+1)*2^(B+1))
// <=> ab ∈ [2^(A+B),2^(A+B+2))
// case 1: A+B <= 62
// 2^(A+B),2^(A+B+2) both increase as A+B increases:
// ab ∈ [2^(A+B),2^(A+B+2))
// ab ∈ [2^(0),2^(62+2))
// <=> ab < 2^64 (note: ab<=MAX_U64: all u64 are < 2^64)
// case 2: A+B >= 64
// ab ∈ [2^64,+∞)
// <=> ab >= 2^64 (note: ab>MAX_U64: no u64 is >= 2^64)
// case 3: A+B = 63
// ab ∈ [2^63,2^65)
// compute the carry: at this point the carry can only be 1 or 0 (bounded by ab<2^65)
// if carry==1
// then ab >= 2^64
// computing carry using (highA = a>>1) instead of (highA = a>>32) allows us some optimizations:
// ((a>>1) * b) doesn't overflow:
// ab ∈ [2^63,2^65)
// (a>>1) * b ∈ [2^62,2^64)
// carry:
// (((a&1)*(b&1)) + (((a>>1)*(b&1))<<1) + (((b>>1)*(a&1))<<1) + (((a>>1)*(b>>1))<<2)) >> 64
// <=> ((((a&1)*(b&1))>>1) + ((a>>1)*(b&1)) + ((b>>1)*(a&1)) + (((a>>1)*(b>>1))<<1)) >> 63
// since: (a>>1)*b == (((a>>1)*(b>>1))<<1) + ((a>>1)*(b&1))
// <=> (((a>>1)*b) + (((a&1)*(b&1))>>1) + ((b>>1)*(a&1))) >> 63
// since: (((a&1)*(b&1))>>1) == 0 : [0,1] * [0,1] ∈ [0,1]; [0,1] >> 1 ∈ [0,0]
// <=> (((a>>1)*b) + ((b>>1)*(a&1))) >> 63
// if we assume the order of likelihood as: A+B < 63, A+B > 63, A+B == 63
// then the branching looks like:
// overflows = A+B >= 63 && (A+B > 63 || (((a>>1)*b) + ((b>>1)*(a&1))))
// since: log2(a) = 63 - clz(a)
// A+B = 126 - clz(a) - clz(b)
// A+B >= 63
// becomes
// 126 - clz(a) - clz(b) >= 63
// -126 + clz(a) + clz(b) <= -63
// clz(a) + clz(b) <= 63
// A+B > 63
// becomes
// clz(a) + clz(b) < 63
int Y = __builtin_clzll(a) + __builtin_clzll(b);
return Y <= 63 && (Y < 63 || ((((a >> 1)*b) + ((b>>1)*(a&1))) >> 63));
}
Если вы просто хотите обнаружить переполнение, как насчет преобразования в удвоение, выполнения умножения и если
| Х | < 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;
}