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


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

#include "TestTime.h"
TestTime tt;
int main(int argc, char* argv[])
{
//char *str=new char[200001];
//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;
}
};