我应该用什么算法来进行高性能的大整数除法?
What algorithm should I use for high-performance large integer division?
我正在将大整数编码为size_t
数组。我已经有了其他的运算(加,减,乘);也可以除以一个个位数。但如果可能的话,我想匹配我的乘法算法的时间复杂度(目前是Toom-Cook)。
我收集了线性时间算法来取我的红利的乘法逆的各种概念。这意味着理论上我可以用与乘法相同的时间复杂度实现除法,因为线性时间操作无论如何通过比较都是"微不足道的"。
我的问题是,我该怎么做呢?哪种类型的乘法逆在实践中是最好的?对64^digitcount
取模?当我用除数乘以乘法逆时,我能避免计算由于整数截断而被丢弃的那部分数据吗?谁能提供C或c++伪代码或给出一个精确的解释应该如何做到这一点?
或者是否存在比基于逆的方法更好的专用除法算法?
编辑:我挖掘了上面提到的"反向"方法。在"计算机编程艺术,卷2:半数值算法"的312页上,Knuth提供了"算法R",这是一个高精度的倒数。他说它的时间复杂度比乘法要小。然而,将其转换为C并对其进行测试是非常重要的,并且在我编写代码之前不清楚将消耗多少开销内存等,这将花费一段时间。如果没人抢在我前面,我就把它贴出来。
对于好的算法,GMP库通常是一个很好的参考。他们记录的除法算法主要依赖于选择一个非常大的基数,这样你就可以用一个4位数除以一个2位数,然后进行长除法。
长除法将需要计算2位数除以1位数的商;这既可以递归地完成,也可以通过预先计算逆并估计商,就像使用Barrett约简一样。
当2n
位数除以n
位数时,递归版本的代价为O(M(n) log(n))
,其中M(n)
为n
位数相乘的代价。
如果你使用牛顿算法计算逆,使用巴雷特约简的版本将花费O(M(n))
,但根据GMP的文档,隐藏常数要大得多,所以这种方法只适用于非常大的除法。
更详细地说,大多数除法算法背后的核心算法是一个"估计商与约简"计算,计算(q,r)
,使
x = qy + r
,但不受0 <= r < y
的限制。典型的循环是
- 估算
x/y
的商q
- 计算相应的减排量
r = x - qy
- 可选地调整商,使减少
r
在某个期望的间隔 - 如果
r
太大,然后用r
代替x
。
x/y
的商将是所有产生的q
的和,r
的最终值将是真正的余数。
分治法通过计算x'/y'
来估计x/y
的商,其中x'
和y'
是x
和y
的前导数字。通过调整它们的大小有很大的优化空间,但是如果x'
的位数是y'
的两倍,那么IIRC可以得到最好的结果。
在我看来,如果坚持整数运算,乘逆方法是最简单的。基本方法是
- 用
- 用
q = 2^(i+j-k) floor(floor(x / 2^i) m / 2^j)
估算x/y
m = floor(2^k / y)
估算y
的反比事实上,实际实现可以容忍m
中的额外错误,如果这意味着你可以使用更快的互惠实现。
错误是一个痛苦的分析,但如果我回忆一下做它的方法,你想选择i
和j
,以便x ~ 2^(i+j)
由于错误是如何积累的,你想选择x / 2^i ~ m^2
以最小化整体工作。
随后的缩减将有r ~ max(x/m, y)
,因此这给出了选择k
的经验法则:您希望m
的大小与每次迭代计算的商的位数有关—或者等效地,每次迭代您想从x
中删除的位数。
我不知道乘法逆算法,但听起来像是蒙哥马利还原法或巴雷特还原法的修改。
我做bigint除法有点不同
参见bignum除法。特别要看一下近似除法和这里的两个链接。一个是我的定点除法,另一个是具有测量的快速乘法算法(如karatsuba,Schönhage-Strassen在NTT上),以及我非常快速的32位基础NTT实现的链接。
我不确定逆乘法是否正确。
主要用于除法器为常数的模运算。我担心对于任意除法,获得bigint逆所需的时间和操作可能比标准除法本身更大,但由于我不熟悉它,我可能错了。
我在实现中看到的最常用的除法是牛顿-拉夫森除法,它与上面链接中的近似除法非常相似。
近似/迭代除法器通常使用乘法来定义其速度。
对于足够小的数字,通常是长二进制除法和32/64位十进制除法,如果不是最快的话,速度足够快:通常它们的开销很小,并且让n
作为处理的最大值(而不是位数!)
二进制除法示例:
是O(log32(n).log2(n)) = O(log^2(n))
.
它循环遍历所有有效位。在每次迭代中,您需要compare, sub, add, bitshift
。这些操作都可以在log32(n)
中完成,log2(n)
是比特数。
下面是我的一个bigint模板(c++)的二进制除法示例:
template <DWORD N> void uint<N>::div(uint &c,uint &d,uint a,uint b)
{
int i,j,sh;
sh=0; c=DWORD(0); d=1;
sh=a.bits()-b.bits();
if (sh<0) sh=0; else { b<<=sh; d<<=sh; }
for (;;)
{
j=geq(a,b);
if (j)
{
c+=d;
sub(a,a,b);
if (j==2) break;
}
if (!sh) break;
b>>=1; d>>=1; sh--;
}
d=a;
}
N
是用于存储bigint型数的32位DWORD
的个数。
-
c = a / b
-
d = a % b
-
qeq(a,b)
是比较:a >= b
大于等于(在log32(n)=N
中完成)
它返回a < b
的0
,a > b
的1
,a == b
的2
-
sub(c,a,b)
是c = a - b
速度提升是通过不使用乘法(如果不计算位移位)获得的
如果你使用像2^32这样的大基数的数字(ALU块),那么你可以在ALU操作中使用32位构建的多项式样式重写整个。
这通常比二进制长除法更快,其思想是将每个DWORD作为单个数字处理,或者递归地将使用的算法除以一半,直到达到CPU能力。参见按半位宽算法除法
最重要的是,当用双位数计算
如果你已经优化了基本操作,那么复杂度可以进一步降低,因为子结果随着迭代变得更小(改变基本操作的复杂度)。一个很好的例子是基于NTT的乘法。
开销会把事情搞砸。
由于这个原因,运行时有时不会复制大O复杂度,因此您应该始终测量阈值并使用更快的方法来使用位计数以获得最大性能并优化您所能。