我正在寻找一种有效(可选标准、优雅且易于实现)的解决方案来乘以相对较大的数字,并将结果存储到一个或多个整数中:
假设我有两个 64 位整数声明如下:
uint64_t a = xxx, b = yyy;
当我这样做时a * b
,如何检测操作是否导致溢出并在这种情况下将进位存储在某处?
请注意,我不想使用任何大型库,因为我对存储数字的方式有限制。
我正在寻找一种有效(可选标准、优雅且易于实现)的解决方案来乘以相对较大的数字,并将结果存储到一个或多个整数中:
假设我有两个 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 ((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" else "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 |
+------------------+----------+----------+----------+----------+----------+
基于这些结果,我们可以得出几个结论:
当 a == 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...
}
问题不仅仅是部分产品,而是任何总和都可能溢出的事实。
如果我必须真正做到这一点,我会用本地汇编语言编写一个扩展乘法例程。 也就是说,例如,将两个 64 位整数相乘得到一个 128 位的结果,该结果存储在两个 64 位寄存器中。所有合理的硬件都在单个本机乘法指令中提供此功能——它不仅可以从 C 中访问。
这是最优雅和最容易编程的解决方案实际上是使用汇编语言的少数情况之一。但它肯定不是便携的:-(
或许解决这个问题最好的方法是有一个函数,将两个 UInt64 相乘,得到一对 UInt64,一个上部和下部的 UInt128 结果。这是解决方案,包括一个以十六进制显示结果的函数。我猜你可能更喜欢 C++ 解决方案,但我有一个有效的 Swift-Solution 显示如何管理问题:
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))
GNU 可移植性库 (Gnulib) 包含一个模块intprops,它具有可有效测试算术运算是否会溢出的宏。
例如,如果发生乘法溢出,INT_MULTIPLY_OVERFLOW (a, b)
将产生1
.
这几天我一直在处理这个问题,我不得不说它给我留下了深刻的印象,我看到人们说知道是否存在溢出的最好方法是划分结果,这完全是低效的,而且不必要。这个函数的要点是它必须尽可能快。
溢出检测有两种选择:
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º- 不可能创建两倍于乘数变量的结果变量:那么您应该使用 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);
}
}
我希望这段代码可以帮助您拥有一个非常高效的程序,并且我希望代码清晰,如果不是,我会发表一些评论。
最好的祝福。
有一个尚未提及的简单(通常非常快速的解决方案)。该解决方案基于这样一个事实,即 n 位乘以 m 位乘法对于 n+m 位或更高的乘积宽度永远不会溢出,但对于小于 n+m-1 的所有结果宽度都会溢出。
因为我的旧描述可能对某些人来说太难阅读了,所以我再试一次:你需要检查两个操作数的前导零的总和。在数学上很容易证明。设 x 为 n 位,y 为 m 位。z = x * y
是 k 位。因为乘积最多可以有 n+m 位大,所以它可能会溢出。比方说。x*y
是 p 位长(没有前导零)。产品的前导零是clz(x * y) = n+m - p
。clz 的行为类似于 log,因此:
clz(x * y) = clz(x) + clz(y) + c with c = either 1 or 0
.
(感谢您在评论中的 c = 1 建议!)
它在k < p <= n+m <=> n+m - k > n+m - p = clz(x * y)
.
现在我们可以使用这个算法:
if max(clz(x * y)) = clz(x) + clz(y) +1 < (n+m - k) --> overflow
if max(clz(x * y)) = clz(x) + clz(y) +1 == (n+m - k) --> overflow if c = 0
else --> no overflow
如何在中间情况下检查溢出?我假设,你有一个乘法指令。然后我们可以很容易地使用它来查看结果的前导零,即:
if clz(x * y / 2) == (n+m - k) <=> msb(x * y/2) == 1 --> overflow
else --> no overflow
您通过将 x/2 视为固定点并将 y 视为普通整数来进行乘法:
msb(x * y/2) = msb(floor(x * y / 2))
floor(x * y/2) = floor(x/2) * y + (lsb(x) * floor(y/2)) = (x >> 1)*y + (x & 1)*(y >> 1)
(这个结果永远不会溢出,以防万一clz(x)+clz(y)+1 == (n+m -k))
诀窍是使用内置/内在函数。在 GCC 中,它看起来是这样的:
static inline int clz(int a) {
if (a == 0) return 32; //only needed for x86 architecture
return __builtin_clz(a);
}
/**@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 = clz(f1) + 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
}
...
只是 n = m = k = 32 位无符号无符号乘法的示例。您可以将其概括为有符号无符号或有符号符号乘法。甚至不需要多位移位(因为一些微控制器仅实现位移位,但有时支持使用一条指令如 Atmega 的乘积除以 2 !)。但是,如果不存在 count-leading-zeroes 指令但长乘法,这可能不会更好。
其他编译器有自己的方式来指定 CLZ 操作的内在函数。与检查乘法的上半部分相比,clz 方法应该比使用高度优化的 128 位乘法检查 64 位溢出更好(在最坏的情况下)。乘法需要线性开销,而计数位只需要线性开销。尝试时,此代码对我来说是开箱即用的。
这是检测两个无符号整数的乘法是否溢出的技巧。
我们观察到,如果我们将 N 位宽的二进制数与 M 位宽的二进制数相乘,则乘积不会超过 N + M 位。
例如,如果要求我们将 3 位数乘以 29 位数,我们知道这不会溢出 32 位。
#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 的内置命令,或者可能与 CPU 的内联汇编一起使用。
现在使用 nsb(a * b) <= nsb(a) + nsb(b) 的事实来检查溢出。当较小时,它总是小 1。
参考 GCC:内置函数:int __builtin_clz (unsigned int x) 返回 x 中前导 0 位的数量,从最高有效位位置开始。如果 x 为 0,则结果未定义。
我今天在想这个,偶然发现了这个问题,我的想法导致了我这个结果。TLDR,虽然我觉得它“优雅”,因为它只使用几行代码(很容易成为单行代码),并且有一些温和的数学,在概念上简化为相对简单的东西,这主要是“有趣的”,我没有没有测试它。
如果您将无符号整数视为基数为 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
这是有道理的。一个数字(例如基数为 10 的 15)环绕的次数是15 / 10 = 1.5
,或 1 次半。但是,我们不能在这里使用 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,则我们不会发生溢出,例如 if 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 时,当我们根据编译器/arch 用于 long double 的有效数字的数量来计算变量时,某些数字可能会出现舍入错误wrap
,我相信应该再多出 20 位才能确定。此结果的另一个问题是,仅通过使用浮点和除法,它的性能可能不如其他一些解决方案。
如果你只是想检测溢出,如何转换为双精度,做乘法,如果
|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;
}