[转贴]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，X(s)=DFT，Y(s)=DFT，那么有，H(s)=X(s)Y(s)
在频域的序列只需要做乘积，就是对应项乘对应项，是O(n)复杂度，因为FFT是一个DFT的快速算法，复杂度在O(nlogn)，所以可以得到同样复杂度的FFT乘法，形式表示成：
h(t)=IFFT[ FFTFFT ]
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
代码现在还在改进加测试（一砣砣的注释代码不敢轻易删），暂不提供。用于测试的代码是：
#include "TestTime.h"
TestTime tt;
int main(int argc, char* argv[])
{
//char *str=new char;
//unsigned long n;
hreal a(0);
double t1k,t10k,t20k,t40k,t80k;
//printf("n=");
//scanf("%u",&n);
t1k=tt.currTime();
a.fibonacci(1000);
t1k=tt.currTime()-t1k;
t10k=tt.currTime();
a.fibonacci(10000);
t10k=tt.currTime()-t10k;
t20k=tt.currTime();
a.fibonacci(20000);
t20k=tt.currTime()-t20k;
t40k=tt.currTime();
a.fibonacci(40000);
t40k=tt.currTime()-t40k;
t80k=tt.currTime();
a.fibonacci(80000);
t80k=tt.currTime()-t80k;
printf("Fibonacci Function
Test\n_______________________\nf(n)\tcost time
(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);
//a.display(str);
//printf("%s\n",str);
//printf("%u! finished.\ncost %lf(ms)\n",n,(t2-t1)*1000);
//delete[] str;
//system("pause");
return 0;
}
TestTime类是用于测时间的，简单一个封装：
#include <windows.h>
class TestTime
{
double freq;
public:
TestTime()
{
LARGE_INTEGER li_freq;
if(!QueryPerformanceFrequency(&li_freq))
freq=-1;
else
{
freq=li_freq.HighPart * 4294967296.0 + li_freq.LowPart;
}
}
double currTime()
{
LARGE_INTEGER nowCount;
QueryPerformanceCounter(&nowCount);
double time=nowCount.HighPart * 4294967296.0 + nowCount.LowPart;
return time/freq;
}
};
来自算法驿站：http://blog.programfan.com/blog.asp?blogid=707 目前实现了，加，减，乘，乘方，阶乘和计算菲数，速度没有大幅度提高，唉，也没那么多时间做了，先放一放。（除法做得不好，太慢了，慢到不能用）最后一次提升的快一些，因为我想到两个实数序列可以合并在一起，做为一个复数的实部和虚部，这样一起做FFT会快一些。推导很简单，公式是：
x(t),y(t)为N点实序列，设X(t)=DFT ,Y(t)=DFT，记w(t)=x(t)+iy(t)，W(t)=DFT，于是由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，开源。
页:
[1]