[转贴]FFT乘法和高精运算库

风花雪月

不准备多写，这些东西网上可以搜出一箩筐来。
高精运算库的构想，是可以快速的应付巨数的运算，包括加，减，乘，除等。比如现在的int是16字节，long32字节，LONGLONG64字节，再大一些就要扩字了，于是就不能用原来汇编给我们用的指令了，需要重新定义运算。
加减很容易，小学加减法，跟手工算一样的，复杂度是O(n)，n表示数据的长度。
乘法如果直接像手工算一样做，复杂度是O(n*n)的，FFT乘法是将数看成多项式：
P1(x)=a0+a1x+a2x^2+...+anx^n，如果取ak在0到9之间，且x取10，那它就表示一个十进制数，如果再用另外一个多项式P2(x)与它相乘，那结果和它所表示的数相乘是一样的数值。
多项式乘法其实质是做两个系数序列的卷积，而FFT正是对这样的卷积做从时域到频域的变化，FFT的本质也是DFT，离散富里叶变换，因为有卷积定理：
如果h(t)=x(t)*(y(t)，且H(s)=DFT[h(t)]，X(s)=DFT[x(t)]，Y(s)=DFT[y(t)]，那么有，H(s)=X(s)Y(s)
在频域的序列只需要做乘积，就是对应项乘对应项，是O(n)复杂度，因为FFT是一个DFT的快速算法，复杂度在O(nlogn)，所以可以得到同样复杂度的FFT乘法，形式表示成：
h(t)=IFFT[ FFT[x(t)]FFT[y(t)] ]
IFFT是快速富里叶反变换，可以转换成FFT，所以一个乘法的时间是3个FFT加一个线性乘积运算，因此也是O(nlogn)。
开始动手做FFT乘法程序之前，先实现一个FFT算法模块，然后要考虑的事有：1，复数如何表示？用float还是double型存实型数？2，进行运算的时候选多大的基？
因为DFT的运算关系到复数，复数好办，用两个实数表示就行了，实数的选取，精度够不够，如何表示精度，误差怎样控制？
如果基是B，参与运算的数最长有N，则一种最极端的情况是，全部取到B-1（测试的时候如果基10，就全9，基100就全99进行测试），做线性卷积回来最大的可能结果是N(B-1)^2，这里要考虑浮点数的精度和需要达到的巨数的长度，很矛盾的选择，一方面想基越大越好，这样大数存起来比较节省，但是基一选大了，精度就会下降，可以实验一下。我实验的结果是float型绝对不够，可能算不是很大的数都错几位，double可以试，但当选B=10000时，N不太大也可能出错，最终选B=100，存储的时候用一个字节存0~99（浪费了一半多，但输出快些）
FFT是技术性最高的，速度也相差很大，有基2的，基4的，混合基的，如果你的计算以乘法为主（比如算阶乘），那FFT会被非常频繁的调用，在这里做一点小小的改动，都会带来速度上的飞跃，可能比你将另外一个地方花几天时间改成汇编要来的经济得多（速度的提高除以你的汗水）。
有了快速乘法，就来实现快速乘幂，比如a^b，通常a是个大数，而b是个不太大的数，算法是将b进行二进制分解，把乘法重叠起来，经过至多2log2(b)次的乘法完成，更准确地说是log2(b)次平方运算+至多log2(b)次乘法运算。你要说平方和乘法不都一样吗，No，平方要快些，a*a，只需要将a变换到频域，自乘，再变回来，2次FFT，而乘法是3次，快一半左右。
除法怎么做？用乘法做，最简单的是用牛顿迭代，设x=1/a，a已知，求x，方程变下形，成1/x-a，导数是-1/x^2，迭代式是x*=x+x-ax^2=x(2-ax)，收敛半径在(0,x*)之间，画图可以观察出来。牛顿法是2阶收敛的，这是说如果收敛的话，每一次迭代的局部误差以双倍的速度减小，比如如果xk的误差是10^-100，那下一次xk+1的误差会成10^-200的数量级，翻倍的收敛。一次迭代要做2次乘法，一次乘法3次FFT，长为n的数据，迭代的次数在log2(n)数量级，那除法的复杂度在O(n(logn)^2)，不太乐观，但比试除的O(n*n)要好。
其实实现了乘法就差不多把整个高精做得差不多了，用你已经实现的再去实现另外其它的东西。你也可以先不去做除法，比如先做开方，或许有更好的方法回来再实现除法。
阶乘是几乎只依赖乘法的模块，做阶乘也挺有意思。阶乘从形式上看是连乘：
n!=1*2*...*n
但是在着手进行连乘之间，希望做什么？使连乘的数变少一些，又想少做几次乘法，又想结果不出错？（鱼和熊掌？）
曾经一道ACM题是算阶乘末尾的0有多少，那个算法会启发你如何从阶乘中提取出素因子的幂，这样一来连乘的项数就会变少，但是遗憾的是如果n比较大素数也相当多，数论里一个经典的完美的证明，就是说这个世界素数是取之不尽数之不完。。。这里还牵涉到计算素数，因为阶乘增涨得比指数还快，数学分析里会让你算一个无穷小，分子是a^n，分母是n!，所以玩阶乘别玩上火，它发起脾气来上帝都受不了（不是只有上帝能追上指数吗），所以通常n较小，可n!很大，在2到n之间求素数，最方便的就是筛法，就是让所有正整数排个队，小的在前面，大的在后面，你从队伍里抓个最小的出来，很容易排头第一个就是，然后宣布他是素数，然后把所有是他的倍数的人踢出队伍，继续筛选，就得到全部的素数，当然当然，1是不算在内的，可以想像你就是1。
阶乘变换成素因子的乘幂，再连乘，就是目前的1.0版，速度测试如下：

Factorial Function Test (1.0)
_______________________
n!      cost time (s)
1k!     0.028163
10k!    0.856568
20k!    2.215239
40k!    6.278020
80k!    16.574058
Test Date\2007\04\11

最后一个实现fibonacci数列地计算，可行的方案有：1，先做开方，用通项求，2，用矩阵运算配合乘法和加法算，我用的后一种，速度比阶乘快多了，结果也没阶乘吓人。

Fibonacci Function Test
_______________________
f(n)    cost time (s)
f(1k)   0.004677
f(10k)  0.058493
f(20k)  0.150679
f(40k)  0.286142
f(80k)  0.749612

代码现在还在改进加测试（一砣砣的注释代码不敢轻易删），暂不提供。用于测试的代码是：

1. #include "TestTime.h"
   
2. TestTime tt;
   
3. int main(int argc, char* argv[])
   
4. {
   
5. //char *str=new char[200001];
   
6. //unsigned long n;
   
7. hreal a(0);
   
8. double t1k,t10k,t20k,t40k,t80k;
   
9. //printf("n=");
   
10. //scanf("%u",&n);
    
11. t1k=tt.currTime();
    
12. a.fibonacci(1000);
    
13. t1k=tt.currTime()-t1k;
    
14. t10k=tt.currTime();
    
15. a.fibonacci(10000);
    
16. t10k=tt.currTime()-t10k;
    
17. t20k=tt.currTime();
    
18. a.fibonacci(20000);
    
19. t20k=tt.currTime()-t20k;
    
20. t40k=tt.currTime();
    
21. a.fibonacci(40000);
    
22. t40k=tt.currTime()-t40k;
    
23. t80k=tt.currTime();
    
24. a.fibonacci(80000);
    
25. t80k=tt.currTime()-t80k;
    
26. printf("Fibonacci Function
    
27. Test\n_______________________\nf(n)\tcost time
    
28. (s)\nf(1k)\t%lf\nf(10k)\t%lf\nf(20k)\t%lf\nf(40k)\t%lf\nf(80k)\t%lf\n",t1k,t10k,t20k,t40k,t80k);
    
29. //a.display(str);
    
30. //printf("%s\n",str);
    
31. //printf("%u! finished.\ncost %lf(ms)\n",n,(t2-t1)*1000);
    
32. //delete[] str;
    
33. //system("pause");
    
34. return 0;
    
35. }

复制代码

TestTime类是用于测时间的，简单一个封装：

1. #include <windows.h>
   
2. class TestTime
   
3. {
   
4. double freq;
   
5. public:
   
6. TestTime()
   
7. {
   
8.   LARGE_INTEGER li_freq;
   
9.   if(!QueryPerformanceFrequency(&li_freq))
   
10.    freq=-1;
    
11.   else
    
12.   {
    
13.    freq=li_freq.HighPart * 4294967296.0 + li_freq.LowPart;
    
14.   }
    
15. }
    
16. double currTime()
    
17. {
    
18.   LARGE_INTEGER nowCount;
    
19.   QueryPerformanceCounter(&nowCount);
    
20.   double time=nowCount.HighPart * 4294967296.0 + nowCount.LowPart;
    
21.   return time/freq;
    
22. }
    
23. };

复制代码

来自算法驿站：http://blog.programfan.com/blog.asp?blogid=707

风花雪月

目前实现了，加，减，乘，乘方，阶乘和计算菲数，速度没有大幅度提高，唉，也没那么多时间做了，先放一放。（除法做得不好，太慢了，慢到不能用）最后一次提升的快一些，因为我想到两个实数序列可以合并在一起，做为一个复数的实部和虚部，这样一起做FFT会快一些。推导很简单，公式是：
x(t),y(t)为N点实序列，设X(t)=DFT[x(t)] ,Y(t)=DFT[y(t)]，记w(t)=x(t)+iy(t)，W(t)=DFT[w(t)]，于是由DFT的线性性质，有W(t)=X(t)+iY(t)，但是不一定就是W(t)的实部和虚部，为了由W(t)得到X(t)和Y(t)，是计算X(t)=(W(t)+W*((N-t))N)/2和Y(t)=(W(t)-W*((N-t))N)/2i。
这样一来，做一次乘法就只需要2次fft。
我的机子比较慢，测了几组数据，结果是：

Factorial Function Test (1.0)
_______________________
n!      cost time (s)
1k!     0.028163
10k!    0.856568
20k!    2.215239
40k!    6.278020
80k!    16.574058
Test Date\2007\04\11
Factorial Function Test (1.01) [增加快速地址变换]
_______________________
n!      Cost Time (s)   Approximations
1000    0.025191        4.023872600770937735437 E 2567
10000   0.735035        2.846259680917054518906 E 35659
20000   1.960328        1.819206320230345134827 E 77337
40000   5.820808        2.091692422212132363320 E 166713
80000   15.231613       3.09772225166224928639 E 357506
Test Date\2007\04\12
Factorial Function Test (1.02) [时频结合免去地址变换]
_______________________
n!      Cost Time (s)   Approximations
1000    0.024248        4.023872600770937735437 E 2567
10000   0.724076        2.846259680917054518906 E 35659
20000   1.910845        1.819206320230345134827 E 77337
40000   5.747496        2.091692422212132363320 E 166713
80000   15.012617       3.09772225166224928639 E 357506
Factorial Function Test (1.03) [代码简单优化]
_______________________
n!      Cost Time (s)   Approximations
1000    0.033531        4.023872600770937735437 E 2567
10000   0.726679        2.846259680917054518906 E 35659
20000   1.770394        1.819206320230345134827 E 77337
40000   5.233378        2.091692422212132363320 E 166713
80000   13.751869       3.09772225166224928639 E 357506
Factorial Function Test (1.03)
_______________________
n!      Cost Time (s)   Approximations
1000    0.044938        4.023872600770937735437 E 2567
10000   0.683504        2.846259680917054518906 E 35659
20000   1.796226        1.819206320230345134827 E 77337
40000   5.204204        2.091692422212132363320 E 166713
80000   14.106093       3.09772225166224928639 E 357506
Factorial Function Test (1.04) [结合硬乘法]
_______________________
n!      Cost Time (s)   Approximations
1000    0.020884        4.023872600770937735437 E 2567
10000   0.699821        2.846259680917054518906 E 35659
20000   1.811458        1.819206320230345134827 E 77337
40000   5.331708        2.091692422212132363320 E 166713
80000   14.116834       3.09772225166224928639 E 357506
Factorial Function Test (1.05) [实数FFT]
_______________________
n!      Cost Time (s)   Approximations
1000    0.027728        4.023872600770937735437 E 2567
10000   0.594507        2.846259680917054518906 E 35659
20000   1.434080        1.819206320230345134827 E 77337
40000   4.070405        2.091692422212132363320 E 166713
80000   10.732933       3.09772225166224928639 E 357506
另外还有算菲数的：
Fibonacci Function Test (1.05)
_______________________
f(n)    Cost Time (s)   Approximations
1000    0.002684        4.34665576869374564356 E 208
10000   0.053351        3.364476487643178326662 E 2089
20000   0.088223        2.531162323732361242240 E 4179
40000   0.217730        1.432600165457807343833 E 8359
80000   0.500366        4.58917898454169424284 E 16718

在网上搜到一些牛人做的，还有功能很完善的HugeCalc，那些速度实在是快，比不得。等以后有时间了，再慢慢玩吧。
截图：




下载包：
http://upload.programfan.com/upfile/200704141621349.rar
只有两个文件：一个hugeint.h一个hugeint.cpp，开源。