多精度计算里,多精度乘法是其中最重要的运算之一,编写的多精度库(类)的其中一个重要效率标识就是其乘法的速度.
根据曾经写的大数类的记忆,简单记录下其中的一些技巧,以备查询.
一般的算法有如:
多精度乘法,所知的几个主要的优化方法有:
1:直接乘法.
2:comba乘法
3:Karatsuba乘法
4:toom_cook乘法
5:FFT乘法
6:FNT算法,或其他类似的有限域上对应类似的FFT算法
7:自己未曾理解的算法.
8:改进技巧,
1:直接乘法.
a b c
* e f g
------------------------
ag bg cg
af bf cf
+ ae be ce
------------------------
(我个人通常称它为硬乘法,后面直接使用这名字),这个是小学就已经学到的普通乘法,一般是没有其他更好的优化条件下使用.一般不大于十几或几十之内都差不多是最好的.
其是n^2复杂度的,准确地说,是C(2,N)
2:comba乘法
a b c
* e f g
------------------------
ag bg cg
af bf cf
+ ae be ce
------------------------
cg = C
bg + cf = D
ag + bf + ce = E
af + be = F
ae = G
-------------------------
实际这个算法和硬乘法的基本原理是一致的,也是n^2复杂度,准确地说也是C(2,N)
但是这里有个特殊的地方,这里计算C,D,E,F,G是可以并行的,这个特性在可允许并行计算时是可以更快地处理.而且这个特性在后面高级优化会有特别用处.
3:Karatsuba乘法
R为进制
a b
* c d
--------------
acR^2 - ((a - b) * (c - d) - ac - bd)R + bd
这里只需要计算 ac,(a - b) *(c - d),bd 3项,较硬乘法里面的4项,少,
其复杂度是n^1.58左右
4:toom_cook乘法
toom_cook乘法,实际可以看成Karatsuba的泛化.
其中可以细化可以有toom3,toom4.toom5.toom6.....
a b c
* d e f
----------
这里abc,def可以看成多项式F(x),G(x)的系数
ax^2 + bx + c
dx^2 + ex + f
-----------------
Ax^4 + Bx^3 + Cx^2 + Dx + E = Q(x)
这里分别取
x = 0,-2,1,-1,2,(或者把2或-2换成无穷大)5个点
这里5个点Q(x)刚好够组成5个等式(也可以看成乘于一个矩阵)
(矩阵Q)|0 0 0 0 1| |F(0)*G(0)|
|16 8 4 2 1| |F(2)*G(2)|
(A,B,C,D,E)* |1 1 1 1 1| = |F(1)*G(1)|
|1 -1 1 -1 1| |F(-1)*G(-1)|
|16 -8 4 -2 1| |F(2)*G(2)|
通过这方程组就能解出A,B,C,D,E,也可以看成求这矩阵的逆.同样可得A,B,C,D,E
toom4,5,6...这些都基本类似方式,选取不同点而得到结果.
至于选取什么点,只是计算技巧问题.
其复杂度是C n^k, k = log(5) / log(3)
其他的 k = log(2k - 1) / log 3;
其中准确的递推公式为 T(3N) = 5T(N) + pN,
由于里面还需要包含了大约5 - 7次左右的N*1位的乘除法法,和若干次N位的加减法,因此此算法额外开销比Karatsuba算法要大,因而尽管复杂度比Karatsuba低,但也只能应用于比Karatsuba更大一些范围才可能有比较好的速度处理.
至于对x选点的问题,一般是看矩阵Q的逆如何更简单一些(不是说数据足够小).
5:FFT乘法
对TOOM_COOK算法选点选择更特殊的点比如是2次幂的根,那实际矩阵Q就变成了范德蒙德矩阵,
这样就利用FFT变换来处理.FFT变换原理不在这展开.
其复杂度是nlogn,对于基2的FFT来说,其比较准确的复杂度是3nlogn 次复数乘法,
6:FNT算法
这里一般是指快速数论变换,而不是快速数论变换里的费马变换.
其复杂度是nlogn,对于基2的FNT来说,其比较准确的复杂度是3 nlogn 次模M剩余系下的乘法运算.
当然还有如二次域,分圆域,伽罗瓦域等一些有限域的扩张下都有类似的处理.通过额外的操作可能换取得到对根和模M的选取更为宽松,
7:自己未曾理解的算法.
自1971年至2007年，最快的是Schnhage–Strassen algorithm，简称SSA,
2007：Fürer‘s algorithm，比SSA更快
这两个似乎都是在有限域上的处理,可惜没看懂.
8:改进技巧,
一般来说
FFT需要3nlogn 次复数乘法,实际应该是3nlogn * 4 = 12nlogn次double乘法
FNT需要3 nlogn 次模M运算下的乘法运算,如果选择一个较大的M 保证计算不溢出,因此M需要选择至少是接近80+ --- 96bit位即3个32bit长度,才能保证计算不出错.即每次运算都需要计算一个96 * 96bit乘法
和一个192 和 96 bit 的除法,才能计算出一次计算结果.实际这样的处理并不是很方便,甚至可能还不如直接计算12 nlogn次double乘法.
因此一般的利用CRT中国剩余定理,
选择合适的 M = M1 * M2 * M3, 其中M1,M2,M3都是32bit范围下的数,这样就变成了
3 nlogn次 模M下的乘法运算,转换成 3 * 3 nlogn = 9 nlogn 次32bit范围内的乘法运算.
CRT最后的合并工作是 n线性的,其工作量相对来nlogn说基本可以忽略
这样处理的优点:
这个较FFT系数上就已经减少了.而且还不算double运算和unsigned int在模乘法速度的差异,
还有一个重要的原因.
FNT并行性要比FFT要好.
理论上
FFT,能够并行处理的是前面两次的变换可以并行处理,
并行处理后
时间上相当于是2 * nlogn次复数乘法 大约是 2 * nlogn * 4 = 8 nlogn次double乘法.
FNT如果能够完全能够做做到并行,前面6次变换是可以并行处理的,后面三次也是完全可以做到并行处理的.即只需要2 nlogn 次乘法,大约是 2nlogn 次unsigned int 模乘法.
当然这个是理想状态,实际虽然有速度差异,但一般没有这么大的区别.