algorithm - 哪个是多精度乘法最有效的算法？

Question

我正在研究一个本地 C++/CLI 类，它执行具有多精度值的整数运算。单个整数由 64 位无符号整数数组表示。符号由 bool 值表示，负值与其绝对值一起存储，而不是作为二进制补码。这使得处理符号问题变得更加容易。目前我正在优化乘法运算。我已经进行了几轮优化，但我的函数仍然需要 的时间是两个 .NET BigInteger 值 的 * 运算符的两倍，这表明进一步优化仍有相当大的潜力。

在寻求帮助之前，让我向您展示我已经尝试过的内容。我的第一次尝试是一种幼稚的方法:使用基本的 64 到 128 位乘法将所有 64 位项目对相乘，然后移位/相加结果。我没有在这里显示代码，因为它非常慢。下一次尝试是 递归分治算法 ，结果证明要好得多。在我的实现中，两个操作数都在中间递归拆分，直到剩下两个 64 位值。这些相乘产生一个 128 位的结果。收集到的基本结果在递归层上一直移位/相加，以产生最终结果。该算法可能受益于需要计算的 64 到 128 位元积少得多的事实，这似乎是主要瓶颈。

所以这是我的代码。第一个片段显示了顶级入口点:

// ----------------------------------------------------------------------------
// Multi-precision multiplication, using a recursive divide-and-conquer plan:
// Left  split: (a*2^k + b)i = ai*2^k + bi
// Right split: a(i*2^k + j) = ai*2^k + aj

public: static UINT64* Mul (UINT64* pu8Factor1,
                            UINT64* pu8Factor2,
                            UINT64  u8Length1,
                            UINT64  u8Length2,
                            UINT64& u8Product)
    {
    UINT64* pu8Product;

    if ((u8Length1 > 0) && (u8Length2 > 0))
        {
        pu8Product = _SnlMemory::Unsigned ((u8Length1 * u8Length2) << 1);

        u8Product  = Mul (pu8Product, pu8Factor1, 0, u8Length1,
                                      pu8Factor2, 0, u8Length2);
        }
    else
        {
        pu8Product = _SnlMemory::Unsigned (0);
        u8Product  = 0;
        }
    return pu8Product;
    }

这些因子作为 UINT64* 数组指针传入，长度分别指定为相应数组中 UINT64 项的数量。该函数分配了一个足够大的内存块来保存最大预期长度的值，它也用作临时从属结果的暂存器。该函数调用另一个 Mul 函数执行递归评估并返回最终结果实际使用的 UINT64 项的数量。

这是分治算法的递归“除法”部分:

// ----------------------------------------------------------------------------
// Recursively expand the arbitrary-precision multiplication to the sum of a
// series of elementary 64-to-128-bit multiplications.

private: static UINT64 Mul (UINT64* pu8Product,
                            UINT64* pu8Factor1,
                            UINT64  u8Offset1,
                            UINT64  u8Length1,
                            UINT64* pu8Factor2,
                            UINT64  u8Offset2,
                            UINT64  u8Length2)
    {
    UINT64 *pu8Lower, u8Lower, *pu8Upper, u8Upper, u8Split;
    UINT64 u8Product = 0;

    if (u8Length1 > 1)
        {
        // left split: (a*2^k + b)i = ai*2^k + bi
        u8Split = u8Length1 >> 1;

        u8Lower = Mul (pu8Lower = pu8Product,
                       pu8Factor1, u8Offset1, u8Split,  // bi
                       pu8Factor2, u8Offset2, u8Length2);

        u8Upper = Mul (pu8Upper = pu8Product + ((u8Split * u8Length2) << 1),
                       pu8Factor1, u8Offset1 + u8Split, // ai
                                   u8Length1 - u8Split,
                       pu8Factor2, u8Offset2, u8Length2);

        u8Product = Mul (u8Split, pu8Lower, u8Lower, pu8Upper, u8Upper);
        }
    else if (u8Length2 > 1)
        {
        // right split: a(i*2^k + j) = ai*2^k + aj
        u8Split = u8Length2 >> 1;

        u8Lower = Mul (pu8Lower = pu8Product,
                       pu8Factor1, u8Offset1, u8Length1, // aj
                       pu8Factor2, u8Offset2, u8Split);

        u8Upper = Mul (pu8Upper = pu8Product + ((u8Length1 * u8Split) << 1),
                       pu8Factor1, u8Offset1, u8Length1, // ai
                       pu8Factor2, u8Offset2 + u8Split,
                                   u8Length2 - u8Split);

        u8Product = Mul (u8Split, pu8Lower, u8Lower, pu8Upper, u8Upper);
        }
    else // recursion base: 64-to-128-bit multiplication
        {
        AsmMul1 (pu8Factor1 [u8Offset1],
                 pu8Factor2 [u8Offset2],
                 u8Lower, u8Upper);

        if (u8Upper > 0)
            {
            pu8Product [u8Product++] = u8Lower;
            pu8Product [u8Product++] = u8Upper;
            }
        else if (u8Lower > 0)
            {
            pu8Product [u8Product++] = u8Lower;
            }
        }
    return u8Product;
    }

在第一个条件分支中，左操作数被拆分。在第二个中，右操作数被拆分。第三个分支是递归基，它调用基本乘法例程:

; -----------------------------------------------------------------------------
; 64-bit to 128-bit multiplication, using the x64 MUL instruction

AsmMul1 proc ; ?AsmMul1@@$$FYAX_K0AEA_K1@Z

; ecx  : Factor1
; edx  : Factor2
; [r8] : ProductL
; [r9] : ProductH

mov  rax, rcx            ; rax = Factor1
mul  rdx                 ; rdx:rax = Factor1 * Factor2
mov  qword ptr [r8], rax ; [r8] = ProductL
mov  qword ptr [r9], rdx ; [r9] = ProductH
ret

AsmMul1 endp

这是一个简单的 ASM PROC，它使用 CPU MUL 指令进行 64 到 128 位乘法运算。我在 ASM 和 C++ 中尝试了其他几个候选者，而这个是最有效的一个。

最后一部分是分治算法的“征服”部分:

// ----------------------------------------------------------------------------
// Shift-add recombination of the results of two partial multiplications.

private: static UINT64 Mul (UINT64  u8Split,
                            UINT64* pu8Lower,
                            UINT64  u8Lower,
                            UINT64* pu8Upper,
                            UINT64  u8Upper)
    {
    FLAG   fCarry;
    UINT64 u8Count, u8Lower1, u8Upper1;
    UINT64 u8Product = u8Lower;

    if (u8Upper > 0)
        {
        u8Count = u8Split + u8Upper;
        fCarry  = false;

        for (u8Product = u8Split; u8Product < u8Count; u8Product++)
            {
            u8Lower1 = u8Product < u8Lower ? pu8Lower [u8Product] : 0;
            u8Upper1 = pu8Upper [u8Product - u8Split];

            if (fCarry)
                {
                pu8Lower [u8Product] = u8Lower1 + u8Upper1 + 1;
                fCarry = u8Lower1 >= MAX_UINT64 - u8Upper1;
                }
            else
                {
                pu8Lower [u8Product] = u8Lower1 + u8Upper1;
                fCarry = u8Lower1 > MAX_UINT64 - u8Upper1;
                }
            }
        if (fCarry)
            {
            pu8Lower [u8Product++] = 1;
            }
        }
    return u8Product;
    }

这里添加了两个部分结果，第二个操作数左移了相应递归步骤的“拆分因子”。

我花了相当多的时间优化代码以提高速度，并取得了相当大的成功，但现在我已经到了一个地步，我看不到任何进一步的可能性，除了使用完全不同的算法。但是，由于我不是数字技巧方面的专家，所以我被困在这里。

希望对如何改进这个计算有一些很酷的想法......

EDIT 2019-03-26: 好吧，有时最好不要尝试变得聪明......经过几次额外的优化尝试，其中一些甚至适度成功，我尝试编写一个真正的愚蠢版本的乘法这只是利用了 _umul128 和 _addcarry_u64 编译器内在函数的计算能力。代码非常简单:

public: static UINT64* Mul (UINT64* pu8Factor1,
                            UINT64* pu8Factor2,
                            UINT64  u8Length1,
                            UINT64  u8Length2,
                            UINT64& u8Product)
    {
    u8Product = u8Length1 + u8Length2;

    CHAR    c1Carry1, c1Carry2;
    UINT64  u8Offset, u8Offset1, u8Offset2, u8Item1, u8Item2, u8Lower, u8Upper;
    UINT64* pu8Product = _SnlMemory::Unsigned (u8Product);

    if (u8Product > 0)
        {
        for (u8Offset1 = 0; u8Offset1 < u8Length1; u8Offset1++)
            {
            u8Offset = u8Offset1;
            u8Item1  = pu8Factor1 [u8Offset1];
            u8Item2  = 0;
            c1Carry1 = 0;
            c1Carry2 = 0;

            for (u8Offset2 = 0; u8Offset2 < u8Length2; u8Offset2++)
                {
                u8Lower = _umul128 (u8Item1, pu8Factor2 [u8Offset2], &u8Upper);

                c1Carry1 = _addcarry_u64 (c1Carry1, u8Item2, u8Lower,
                                          &u8Item2);

                c1Carry2 = _addcarry_u64 (c1Carry2, u8Item2,
                                          pu8Product  [u8Offset],
                                          pu8Product + u8Offset);
                u8Item2 = u8Upper;
                u8Offset++;
                }
            if (c1Carry1 != 0)
                {
                c1Carry2 = _addcarry_u64 (c1Carry2, u8Item2 + 1,
                                          pu8Product  [u8Offset],
                                          pu8Product + u8Offset);
                }
            else if (u8Item2 != 0)
                {
                c1Carry2 = _addcarry_u64 (c1Carry2, u8Item2,
                                          pu8Product  [u8Offset],
                                          pu8Product + u8Offset);
                }
            }
        if (pu8Product [u8Product - 1] == 0)
            {
            u8Product--;
            }
        }
    return pu8Product;
    }

它在堆上创建一个结果缓冲区，足够大以容纳产品的最大大小，并在两个嵌套循环中结合使用 _umul128 的两个波纹进位加法流进行基本的 64-zo-128 位 _addcarry_u64 乘法。这个版本的性能是迄今为止我尝试过的最好的!它比等效的 .NET BigInteger 乘法快 10 倍左右，所以最终我实现了 20 倍的加速。

Answer 1

正如我们在 reference source 中看到的那样，.NET 中的 BigInteger 使用相当慢的乘法算法，通常使用 32x32->64 乘法的二次时间算法。但它是用低开销编写的:迭代、很少分配，并且不调用不可内联的 ASM 过程。部分产品会立即添加到结果中，而不是单独实现。

不可内联的 ASM 过程可以替换为 _umul128固有的。手动进位计算(条件 +1 和确定输出进位)可以替换为 _addcarry_u64固有的。

更高级的算法，例如 Karatsuba 乘法和 Toom-Cook 乘法可能是有效的，但当递归一直执行到单肢级别时就不行了——这远远超过了开销超过保存的基本乘法的点。作为一个具体的例子，this implementation Java 的 BigInteger 的 80 个肢体切换到 Karatsuba(2560 位，因为它们使用 32 位肢体)，以及 240 个肢体的 3-way Toom-Cook。考虑到 80 的阈值，只有 64 条肢体，我不会期望那里有太多的 yield ，如果有的话。