在64位中进行组合乘法除法运算的最准确方法是什么?

对于在32位和64位程序(在Visual C ++中)都能工作的64位整数,我能够进行乘法除法运算的最准确方法是什么? (如果溢出,我需要结果mod 2 64.

(我正在寻找类似MulDiv64的东西,除了这个使用内联汇编,它只适用于32位程序。)

显然,可以进行double和后退,但我想知道是否有更准确的方法并不太复杂。 (即我不是在寻找任意精度的算术库!)

由于这是标记为Visual C ++,我将提供一个滥用MSVC特定内在函数的解决方案。

这个例子相当复杂。 它是GMP和java.math.BigInteger用于大分区的相同算法的高度简化版本。

虽然我有一个更简单的算法,但它可能慢了约30倍。

此解决方案具有以下约束/行为:

  • 它需要x64。 它不会在x86上编译。
  • 商不为零。
  • 如果它溢出64位,则商饱和。

请注意,这是针对无符号整数的情况。 为此创建一个包装器以使其适用于已签名的案例是微不足道的。 此示例还应生成正确的截断结果。

此代码未经过完全测试。 但是,它已经通过了我抛出的所有测试用例。
(即使我故意构建以试图打破算法的情况。)

 #include  uint64_t muldiv2(uint64_t a, uint64_t b, uint64_t c){ // Normalize divisor unsigned long shift; _BitScanReverse64(&shift,c); shift = 63 - shift; c <<= shift; // Multiply a = _umul128(a,b,&b); if (((b << shift) >> shift) != b){ cout << "Overflow" << endl; return 0xffffffffffffffff; } b = __shiftleft128(a,b,shift); a <<= shift; uint32_t div; uint32_t q0,q1; uint64_t t0,t1; // 1st Reduction div = (uint32_t)(c >> 32); t0 = b / div; if (t0 > 0xffffffff) t0 = 0xffffffff; q1 = (uint32_t)t0; while (1){ t0 = _umul128(c,(uint64_t)q1 << 32,&t1); if (t1 < b || (t1 == b && t0 <= a)) break; q1--; // cout << "correction 0" << endl; } b -= t1; if (t0 > a) b--; a -= t0; if (b > 0xffffffff){ cout << "Overflow" << endl; return 0xffffffffffffffff; } // 2nd reduction t0 = ((b << 32) | (a >> 32)) / div; if (t0 > 0xffffffff) t0 = 0xffffffff; q0 = (uint32_t)t0; while (1){ t0 = _umul128(c,q0,&t1); if (t1 < b || (t1 == b && t0 <= a)) break; q0--; // cout << "correction 1" << endl; } // // (a - t0) gives the modulus. // a -= t0; return ((uint64_t)q1 << 32) | q0; } 

请注意,如果您不需要完美截断的结果,则可以完全删除最后一个循环。 如果这样做,答案将不超过正确商数2。

测试用例:

 cout << muldiv2(4984198405165151231,6132198419878046132,9156498145135109843) << endl; cout << muldiv2(11540173641653250113, 10150593219136339683, 13592284235543989460) << endl; cout << muldiv2(449033535071450778, 3155170653582908051, 4945421831474875872) << endl; cout << muldiv2(303601908757, 829267376026, 659820219978) << endl; cout << muldiv2(449033535071450778, 829267376026, 659820219978) << endl; cout << muldiv2(1234568, 829267376026, 1) << endl; cout << muldiv2(6991754535226557229, 7798003721120799096, 4923601287520449332) << endl; cout << muldiv2(9223372036854775808, 2147483648, 18446744073709551615) << endl; cout << muldiv2(9223372032559808512, 9223372036854775807, 9223372036854775807) << endl; cout << muldiv2(9223372032559808512, 9223372036854775807, 12) << endl; cout << muldiv2(18446744073709551615, 18446744073709551615, 9223372036854775808) << endl; 

输出:

 3337967539561099935 8618095846487663363 286482625873293138 381569328444 564348969767547451 1023786965885666768 11073546515850664288 1073741824 9223372032559808512 Overflow 18446744073709551615 Overflow 18446744073709551615 

你只需要64位整数。 有一些冗余操作但允许在调试器中使用10作为基础和步骤。

 uint64_t const base = 1ULL<<32; uint64_t const maxdiv = (base-1)*base + (base-1); uint64_t multdiv(uint64_t a, uint64_t b, uint64_t c) { // First get the easy thing uint64_t res = (a/c) * b + (a%c) * (b/c); a %= c; b %= c; // Are we done? if (a == 0 || b == 0) return res; // Is it easy to compute what remain to be added? if (c < base) return res + (a*b/c); // Now 0 < a < c, 0 < b < c, c >= 1ULL // Normalize uint64_t norm = maxdiv/c; c *= norm; a *= norm; // split into 2 digits uint64_t ah = a / base, al = a % base; uint64_t bh = b / base, bl = b % base; uint64_t ch = c / base, cl = c % base; // compute the product uint64_t p0 = al*bl; uint64_t p1 = p0 / base + al*bh; p0 %= base; uint64_t p2 = p1 / base + ah*bh; p1 = (p1 % base) + ah * bl; p2 += p1 / base; p1 %= base; // p2 holds 2 digits, p1 and p0 one // first digit is easy, not null only in case of overflow uint64_t q2 = p2 / c; p2 = p2 % c; // second digit, estimate uint64_t q1 = p2 / ch; // and now adjust uint64_t rhat = p2 % ch; // the loop can be unrolled, it will be executed at most twice for // even bases -- three times for odd one -- due to the normalisation above while (q1 >= base || (rhat < base && q1*cl > rhat*base+p1)) { q1--; rhat += ch; } // subtract p1 = ((p2 % base) * base + p1) - q1 * cl; p2 = (p2 / base * base + p1 / base) - q1 * ch; p1 = p1 % base + (p2 % base) * base; // now p1 hold 2 digits, p0 one and p2 is to be ignored uint64_t q0 = p1 / ch; rhat = p1 % ch; while (q0 >= base || (rhat < base && q0*cl > rhat*base+p0)) { q0--; rhat += ch; } // we don't need to do the subtraction (needed only to get the remainder, // in which case we have to divide it by norm) return res + q0 + q1 * base; // + q2 *base*base } 

这是一个社区维基的答案,因为它实际上只是一堆指向其他论文/参考的指针(我无法发布相关代码)。

使用铅笔和纸张技术的简单应用,将两个64位整数乘以128位结果非常容易,每个人都在小学就读。

GregS的评论是正确的:Knuth在第4.3.1节“多精度算术/经典算法”(我的副本中的第255-265页)末尾的“计算机编程艺术,第二版,第2卷/数值算法”中涵盖了部门。 这不是一个简单的阅读,至少不是像我这样的人忘记了7年级代数以外的大多数数学。 就在此之前,Knuth也涵盖了事物的乘法方面。

思想的其他一些选择(这些注释用于除法算法,但大多数也讨论乘法):

  • Jack Crenshaw在嵌入式系统编程杂志1997的一系列文章中以更易读的方式介绍了Knuth分区算法(不幸的是,我的笔记没有确切的问题)。 可悲的是,旧ESP问题的文章在网上很难找到。 如果您可以访问大学图书馆,可能会有一些问题或ESP CD-ROM库的副本可供您使用。
  • 微软研究院的Thomas Rodeheffer发表了一篇关于软件整合部门的论文: http ://research.microsoft.com/pubs/70645/tr-2008-141.pdf
  • KarlHasselström关于“大整数快速划分”的论文: http : //www.treskal.com/kalle/exjobb/original-report.pdf
  • Randall Hyde的“汇编语言艺术”(http://webster.cs.ucr.edu/AoA/Windows/HTML/AoATOC.html),特别是第四卷第4.2.5节(扩展精度部门): http:// webster .cs.ucr.edu / AoA / Windows / HTML / AdvancedArithmetica2.html#998729这是Hyde的x86汇编语言的变体,但也有一些伪代码和足够的解释将算法移植到C.它也很慢 – 执行分区逐位…

你不需要任意精度算术。 您只需要128位算术。 即你需要64 * 64 = 128乘法和128/64 = 64除法(具有适当的溢出行为)。 这并不难以手动实现。

那么你可以将你的64位操作数切换成32位块(低和高部分)。 比你想做的对他们的操作。 所有中间结果都将小于64位,因此可以存储在您拥有的数据类型中。

你是否在VC ++中拥有COMP类型(基于x87的64位整数类型)? 我过去常常在Delphi中使用过它,当时我需要64位整数数学。 多年来,它比基于库的64位整数数学更快 – 当然涉及到分割时。

在Delphi 2007(我安装的最新版本 – 32位)中,我会像这样实现MulDiv64:

 function MulDiv64(const a1, a2, a3: int64): int64; var c1: comp absolute a1; c2: comp absolute a2; c3: comp absolute a3; r: comp absolute result; begin r := c1*c2/c3; end; 

(那些奇怪的绝对语句将comp变量叠加在它们的64位整数计数器部分之上。我会使用简单的类型转换,除了Delphi编译器对此感到困惑 – 可能是因为Delphi语言(或者他们现在称之为的任何语言)类型转换(重新解释)和值类型转换之间没有明确的语法区别。)

无论如何,Delphi 2007将以上内容呈现如下:

 0046129C 55 push ebp 0046129D 8BEC mov ebp,esp 0046129F 83C4F8 add esp,-$08 004612A2 DF6D18 fild qword ptr [ebp+$18] 004612A5 DF6D10 fild qword ptr [ebp+$10] 004612A8 DEC9 fmulp st(1) 004612AA DF6D08 fild qword ptr [ebp+$08] 004612AD DEF9 fdivp st(1) 004612AF DF7DF8 fistp qword ptr [ebp-$08] 004612B2 9B wait 004612B3 8B45F8 mov eax,[ebp-$08] 004612B6 8B55FC mov edx,[ebp-$04] 004612B9 59 pop ecx 004612BA 59 pop ecx 004612BB 5D pop ebp 004612BC C21800 ret $0018 

以下语句产生256204778801521550,这似乎是正确的。

 writeln(MulDiv64($aaaaaaaaaaaaaaa, $555555555555555, $1000000000000000)); 

如果你想将它实现为VC ++内联汇编,你可能需要对默认的舍入标志进行一些调整来完成同样的事情,我不知道 – 我没有必要找出 – 但:)

对于64位模式代码,您可以实现64 * 64 = 128乘法,类似于128/64 = 64:64除法的实现。

对于32位代码,它将更复杂,因为没有CPU指令可以在32位模式下对这些长操作数进行乘法或除法,并且您必须将几个较小的乘法组合成一个较大的乘法并重新实现长除法。

您可以使用此答案中的代码作为构建长分区的框架。

当然,如果你的分频器总是小于2 32 (或者更好的是2 16 ),你可以通过链接几个分区(除以被分红的最高64位(或32位))在32位模式下进行更快的除法。 32位(或16位)除数,然后将该除法的余数与被除数的下一个64位(32)进行组合,并除以除数,并继续这样做直到用完整个被除数)。 此外,如果除数很大但可以考虑到足够小的数,则使用此链式除法除以其因子将优于经典的循环解。

这是一个可以使用的近似方法:(全精度,除非> 0x7FFFFFFF或b> 0x7FFFFFFF且c大于a或b)

 constexpr int64_t muldiv(int64_t a, int64_t b, int64_t c, unsigned n = 0) { return (a < 0x7FFFFFFF && b < 0x7FFFFFFF) ? (a * b) / c : (n != 2) ? (c <= a) ? ((a / c) * b + muldiv(b, a % c, c, n + 1)) : muldiv(a, b, c / 2) / 2 : 0; } 

模数用于找出精度损失,然后将其插回到函数中。 这类似于经典的除法算法。

选择2是因为(x % x) % x = x % x

考虑你想要乘以a乘以b然后除以d

 uint64_t LossyMulDiv64(uint64_t a, uint64_t b, uint64_t d) { long double f = long double(b)/d; uint64_t highPart = uint64_t((a & ~0xffffffff) * f + 0.5); uint64_t lowPart = uint64_t((a & 0xffffffff) * f + 0.5); return highPart + lowPart; } 

此代码将a的值拆分为更高和更低的32位部分,然后将32位部分分别乘以b的52位精确比率d ,对部分乘法进行舍入,并将它们加总回整数。 一些精度仍然会丢失,但结果比仅return a * double(b) / d;更精确return a * double(b) / d;

如果您只需要支持Windows 7及更高版本,一个好方法是:

 #include  #include  #pragma comment( lib, "mfplat.lib" ) uint64_t mulDiv64( uint64_t a, uint64_t b, uint64_t c ) { assert( a <= LLONG_MAX && b <= LLONG_MAX && c <= LLONG_MAX ); // https://docs.microsoft.com/en-us/windows/desktop/api/Mfapi/nf-mfapi-mfllmuldiv return (uint64_t)MFllMulDiv( (__int64)a, (__int64)b, (__int64)c, (__int64)c / 2 ); } 

这种方法比其他答案简单得多,它将结果舍入而不是截断,并且适用于所有Windows平台,包括ARM。