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

Question

在 32 位和 64 位程序（在 Visual C++ 中）中，我可以对 64 位整数进行乘除运算的最准确方法是什么？ （如果溢出，我需要结果 mod 2⁶⁴。）

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

显然，转换为double 并返回是可能的，但我想知道是否有一种不太复杂的更准确的方法。 （即我不是在这里寻找任意精度的算术库！）

Answer 1

您可以将 64 位操作数分成 32 位块（低位和高位）。比对它们进行您想要的操作。所有中间结果都将小于 64 位，因此可以存储在您拥有的数据类型中。

Answer 2

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

Answer 3

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

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

您可以使用this answer 中的代码作为构建长除法的框架。

当然，如果您的除法器总是小于 2³²（或者更好的是 2¹⁶），您可以通过链接在 32 位模式下进行更快的除法几个除法（将被除数的最高有效 64（或 32）位除以 32 位（或 16 位）除数，然后将此除法的余数与被除数的下 64 (32) 位相结合，然后除以除数并继续这样做，直到你用完整个股息）。此外，如果除数很大但可以分解为足够小的数字，则使用这种链式除法除以它的因数将比经典循环解决方案更好。

Answer 4

在 VC++ 中是否有 COMP 类型（基于 x87 的 64 位整数类型）可供使用？过去，当我需要 64 位整数数学时，我偶尔会在 Delphi 中使用它。多年来，它比基于库的 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;

（那些奇怪的 absolute 语句将 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++ 内联汇编，您可能需要对默认舍入标志进行一些调整以完成相同的事情，我不知道 - 我不需要找出来——但:)

Answer 5

这是一种您可以使用的近似方法：（全精度，除非 a > 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。

Answer 6

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

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

GregS 的评论是正确的：Knuth 在第 4.3.1 节多精度算术/经典算法（第 255 - 265 页在我的复制）。这不是一本容易读的书，至少对于像我这样忘记了七年级代数以外的大多数数学的人来说不是。就在之前，Knuth 也涵盖了乘法方面。

其他一些想法选项（这些注释是针对除法算法的，但大多数也讨论乘法）：

• Jack Crenshaw 在 1997 年嵌入式系统编程杂志的一系列文章中以更易读的方式介绍了 Knuth 除法算法（不幸的是，我的笔记没有确切的问题）。可悲的是，旧 ESP 问题的文章在网上并不容易找到。如果您可以访问大学图书馆，您可能可以获得一些过刊或 ESP CD-ROM 图书馆的副本。
• Microsoft Research 的 Thomas Rodeheffer 有一篇关于软件整数部的论文：http://research.microsoft.com/pubs/70645/tr-2008-141.pdf
• Karl Hasselströ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 节（扩展精度划分）：@987654323 @ 这是 Hyde 的 x86 汇编语言的变体，但也有一些伪代码和足够的解释将算法移植到 C。它也很慢 - 逐位执行除法......

Answer 7

由于这被标记为 Visual C++，我将给出一个滥用 MSVC 特定内在函数的解决方案。

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

虽然我想到了一个更简单的算法，但它可能慢了大约 30 倍。

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

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

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

这段代码没有经过全面测试。但是，它已经通过了我提交给它的所有测试用例。
（即使是我故意构建的尝试破坏的用例）算法。）

#include <intrin.h>

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

Answer 8

您只需要 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
}

Answer 9

假设您要将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 与 d 的 52 位精确比率，将部分乘法四舍五入，然后将它们加起来返回一个整数。 一些精度仍然会丢失，但结果比简单的 return a * double(b) / d; 更精确。

Answer 10

如果您只需要支持 Windows 7 和更新版本，一个好方法是：

#include <mfapi.h>
#include <assert.h>
#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 );
}

此方法比此处其他答案中的方法简单得多，它对结果进行舍入而不是截断，并且适用于包括 ARM 在内的所有 Windows 平台。

Answer 11

由于这是标记为 Visual C++，您可以使用 newly available intrinsics:

• 多：_mul128()，_umul128()
• 分区：_div128()，_udiv128()

uint64_t muldiv_u64(uint64_t a, uint64_t b, uint64_t c)
{
    uint64_t highProduct;
    uint64_t lowProduct = _umul128(a, b, &highProduct);
    uint64_t remainder;
    return _udiv128(highProduct, lowProduct, c, &remainder);
}

如果您需要签名的 mul-div，则只需使用不带 u 的版本

另见

• How can I descale x by n/d, when x*n overflows?
• How can I multiply and divide 64-bit ints accurately?
• How to multiply a 64 bit integer by a fraction in C++ while minimizing error?
• (a * b) / c MulDiv and dealing with overflow from intermediate multiplication