哪种算法是最有效的多精度乘法算法？

Question

我正在处理一个本地C++/CLI类，它使用多精度值执行整数运算。单个整数由64位无符号整数数组表示.符号用布尔值表示，负值是用它们的绝对值存储的，而不是作为两个的补充。这使得处理符号问题变得容易多了。目前我正在优化乘法运算。我已经做了几次优化循环，但是函数仍然需要是两个.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程序，它使用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;
    }

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

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

希望能对如何改进这个计算提出一些好的想法.

编辑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;
    }

它在堆上创建一个结果缓冲区，足够大到足以容纳产品的最大大小，并在两个嵌套循环中执行基本的64-zo-128位_umul128乘法，以及使用_addcarry_u64进行两个波纹携带加法。到目前为止，这个版本的性能是我尝试过的最好的！它比等效的.NET BigInteger乘法快10倍，所以最后我达到了20倍的速度。

Answer 1

正如我们在参考源中看到的，BigInteger在.NET中使用了一个相当慢的乘法算法，通常使用32x32->64的二次时间算法。但是它是低开销的:迭代，很少分配，没有调用不可否认的ASM过程。部分产品将立即添加到结果中，而不是单独实现。

不可否认的ASM过程可以用umul128内禀来代替.手动进位计算(条件+1和确定输出进位)都可以由_addcarry_u64内蕴代替。

诸如Karatsuba乘法和Toom-Cook乘法等更高级的算法可能是有效的，但当递归一直到单肢级时就不会有效了--这远远超过了开销超过节省的基本乘法的程度。作为一个具体的例子，这一实现 of Java的BigInteger切换到Karatsuba的80个肢体(2560位，因为他们使用32位肢体)，和3路Toom为240个肢体。鉴于这一门槛为80，只有64条腿，我不指望有太多的收获，如果有。

Answer 2

我实现了Karatsuba作为一个递归子程序，它完成了连续的二分运算，而且比标准的长乘法要快得多。停止2:1-分裂，当你达到一个长度，其中有一个硬件复印机，当然！您还可以通过在两个或多个核心机器上进行多读处理来进行并行处理。