哪种是最有效的多精度乘法算法？
Which is the most efficient algorithm for multiple-precision multiplication?
我正在开发本机 C++/CLI class，它执行具有多精度值的整数运算。单个整数由 64 位无符号整数数组表示。符号由布尔值表示，负值与其绝对值一起存储，而不是二进制补码。这使得处理符号问题变得更加容易。目前我正在优化乘法运算。我已经完成了几轮优化，但我的函数仍然需要 两倍于 * 运算符的时间 ，这表明进一步优化的潜力仍然很大。
在寻求帮助之前，让我向您展示一下我已经尝试过的方法。我的第一次尝试是一种天真的方法：使用基本的 64 位到 128 位乘法将所有 64 位项的对相乘，然后 shift/add 结果。我没有在这里显示代码，因为它非常慢。下一次尝试是递归分治算法，结果证明要好得多。在我的实现中，两个操作数在中间递归拆分，直到剩下两个 64 位值。这些相乘产生 128 位结果。收集到的基本结果一直 shift/add 向上递归层产生最终结果。该算法可能受益于需要计算的 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 函数执行递归计算和 returns 最终结果实际使用的 UINT64 项目数。
这是分而治之算法的递归"divide"部分：
// ----------------------------------------------------------------------------
// 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++ 中尝试过其他几个候选者，这个是最有效的。
最后一部分是分治算法的"conquer"部分：
// ----------------------------------------------------------------------------
// 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;
}
这里添加了两个部分结果，第二个操作数左移相应递归步骤的 "split factor"。
我花了相当多的时间来优化代码以提高速度，并取得了一些成功，但现在我已经到了一个地步，除了使用完全不同的算法之外，我看不到任何进一步的可能性。但是，由于我不是数字技巧方面的专家，所以我被困在这里。
希望获得一些关于如何改进此计算的绝妙想法...
编辑 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 倍的加速。
正如我们在reference source中看到的，.NET中的BigInteger使用了相当慢的乘法算法，通常的二次时间算法使用32x32->64次乘法。但它是用低开销编写的：迭代、很少的分配，并且不调用不可内联的 ASM 过程。部分产品会立即添加到结果中，而不是单独具体化。
不可内联的 ASM 过程可以替换为 _umul128 内在过程。手动进位计算（条件 +1 和确定输出进位）可以替换为 _addcarry_u64 内在计算。
诸如 Karatsuba 乘法和 Toom-Cook 乘法之类的高级算法可能是有效的，但当递归一直执行到单肢级别时就不是这样了——这远远超过了开销超过保存的基本级别的点乘法。作为具体示例，Java 的 BigInteger 的 this implementation 切换到 Karatsuba 80 个肢体（2560 位，因为他们使用 32 位肢体），并切换到 3 向 Toom-Cook 240 个肢体。考虑到 80 的阈值，只有 64 个肢体，我不会期望在那里获得太多收益，如果有的话。
我正在开发本机 C++/CLI class，它执行具有多精度值的整数运算。单个整数由 64 位无符号整数数组表示。符号由布尔值表示，负值与其绝对值一起存储，而不是二进制补码。这使得处理符号问题变得更加容易。目前我正在优化乘法运算。我已经完成了几轮优化，但我的函数仍然需要 两倍于 * 运算符的时间 ，这表明进一步优化的潜力仍然很大。
在寻求帮助之前，让我向您展示一下我已经尝试过的方法。我的第一次尝试是一种天真的方法：使用基本的 64 位到 128 位乘法将所有 64 位项的对相乘，然后 shift/add 结果。我没有在这里显示代码，因为它非常慢。下一次尝试是递归分治算法，结果证明要好得多。在我的实现中，两个操作数在中间递归拆分，直到剩下两个 64 位值。这些相乘产生 128 位结果。收集到的基本结果一直 shift/add 向上递归层产生最终结果。该算法可能受益于需要计算的 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 函数执行递归计算和 returns 最终结果实际使用的 UINT64 项目数。
这是分而治之算法的递归"divide"部分：
// ----------------------------------------------------------------------------
// 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++ 中尝试过其他几个候选者，这个是最有效的。
最后一部分是分治算法的"conquer"部分：
// ----------------------------------------------------------------------------
// 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;
}
这里添加了两个部分结果，第二个操作数左移相应递归步骤的 "split factor"。
我花了相当多的时间来优化代码以提高速度，并取得了一些成功，但现在我已经到了一个地步，除了使用完全不同的算法之外，我看不到任何进一步的可能性。但是，由于我不是数字技巧方面的专家，所以我被困在这里。
希望获得一些关于如何改进此计算的绝妙想法...
编辑 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 倍的加速。
正如我们在reference source中看到的，.NET中的BigInteger使用了相当慢的乘法算法，通常的二次时间算法使用32x32->64次乘法。但它是用低开销编写的：迭代、很少的分配，并且不调用不可内联的 ASM 过程。部分产品会立即添加到结果中，而不是单独具体化。
不可内联的 ASM 过程可以替换为 _umul128 内在过程。手动进位计算（条件 +1 和确定输出进位）可以替换为 _addcarry_u64 内在计算。
诸如 Karatsuba 乘法和 Toom-Cook 乘法之类的高级算法可能是有效的，但当递归一直执行到单肢级别时就不是这样了——这远远超过了开销超过保存的基本级别的点乘法。作为具体示例，Java 的 BigInteger 的 this implementation 切换到 Karatsuba 80 个肢体（2560 位，因为他们使用 32 位肢体），并切换到 3 向 Toom-Cook 240 个肢体。考虑到 80 的阈值，只有 64 个肢体，我不会期望在那里获得太多收益，如果有的话。