algorithm - 哪个是多精度乘法最有效的算法？
问题描述
我正在研究一个本地 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 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;
}
这里添加了两个部分结果，第二个操作数左移了相应递归步骤的“拆分因子”。
我花了相当多的时间优化代码以提高速度，并取得了相当大的成功，但现在我已经到了一个地步，我看不到任何进一步的可能性，除了使用完全不同的算法。但是，由于我不是数字技巧方面的专家，所以我被困在这里。
希望对如何改进这个计算有一些很酷的想法......
编辑 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在两个嵌套循环中结合使用_addcarry_u64. 这个版本的性能是迄今为止我尝试过的最好的！它比等效的 .NET BigInteger 乘法快 10 倍左右，所以最终我实现了 20 倍的加速。
解决方案
正如我们在参考资料中看到的那样，.NET 中的 BigInteger 使用相当慢的乘法算法，通常使用 32x32->64 乘法的二次时间算法。但它是用低开销编写的：迭代、很少分配，并且不调用不可内联的 ASM 过程。部分产品会立即添加到结果中，而不是单独实现。
不可内联的 ASM 过程可以替换为_umul128内在函数。手动进位计算（条件+1进位和确定输出进位）可以替换为_addcarry_u64内在进位。
更高级的算法，如 Karatsuba 乘法和 Toom-Cook 乘法可能是有效的，但当递归一直执行到单肢级别时就不行了——这远远超过了开销超过保存的基本乘法的点。作为一个具体的例子， Java 的 BigInteger 的这个实现切换到 Karatsuba 以获取 80 个肢体（2560 位，因为它们使用 32 位肢体），并切换到 3-way Toom-Cook 以获取 240 个肢体。考虑到 80 的阈值，只有 64 条肢体，我不会期望那里有太多的收益，如果有的话。
推荐阅读
- flutter - Flutter Carousel z 顺序重叠
- javascript - 如何大规模创建/高效制作数千条快速 js 路由？
- google-bigquery - 在 bigquery 中使用 bigquery select 语句的 JSON 格式的现有表的架构
- javascript - 使用嵌套连接创建的 prisma 抛出错误
- sql - 在 SQL 中使用 Self Join 获取金额
- ldap - LDAP - 向 LDAP URL 添加过滤器，以便只有属于某个组的用户才能访问应用程序
- android - 在此处为地图对象设置 zindex 地图 android
- node.js - 根据挂钟时间组合视频和音频 .webm 文件
- java - 通过添加 256 将字节数组中的负字节变为正数 - java
- python - 在没有 pip 的情况下安装 Python 包