空间管理您的位置: 51Testing软件测试网 » 努力学习美好将来 » 日志

FFT乘法和高精运算库(1)

上一篇 / 下一篇 2007-04-18 13:46:32 / 个人分类：算法

!WI `J%]5O\]"LE0高精运算库的构想，是可以快速的应付巨数的运算，包括加，减，乘，除等。比如现在的int是16字节，long32字节，LONGLONG64字节，再大一些就要扩字了，于是就不能用原来汇编给我们用的指令了，需要重新定义运算。51Testing软件测试网IgDqPjZ4E

Up"cy{c*q.uIm0加减很容易，小学加减法，跟手工算一样的，复杂度是O(n)，n表示数据的长度。

U2c1e1m'l le}e0乘法如果直接像手工算一样做，复杂度是O(n*n)的，FFT乘法是将数看成多项式：

]1N&H)O8_0P1(x)=a0+a1x+a2x^2+...+anx^n，如果取ak在0到9之间，且x取10，那它就表示一个十进制数，如果再用另外一个多项式P2(x)与它相乘，那结果和它所表示的数相乘是一样的数值。51Testing软件测试网y!G,{,I%W

多项式乘法其实质是做两个系数序列的卷积，而FFT正是对这样的卷积做从时域到频域的变化，FFT的本质也是DFT，离散富里叶变换，因为有卷积定理：

&~7P2f }rEPQ0如果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)51Testing软件测试网#`7]9RNwm

在频域的序列只需要做乘积，就是对应项乘对应项，是O(n)复杂度，因为FFT是一个DFT的快速算法，复杂度在O(nlogn)，所以可以得到同样复杂度的FFT乘法，形式表示成：

h(t)=IFFT[ FFT[x(t)]FFT[y(t)] ]

IFFT是快速富里叶反变换，可以转换成FFT，所以一个乘法的时间是3个FFT加一个线性乘积运算，因此也是O(nlogn)。51Testing软件测试网M v2C7S'Y1z

7I Ve*X)@4K5nU0开始动手做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会被非常频繁的调用，在这里做一点小小的改动，都会带来速度上的飞跃，可能比你将另外一个地方花几天时间改成汇编要来的经济得多（速度的提高除以你的汗水）。51Testing软件测试网Z`V6u{!Lx

有了快速乘法，就来实现快速乘幂，比如a^b，通常a是个大数，而b是个不太大的数，算法是将b进行二进制分解，把乘法重叠起来，经过至多2log2(b)次的乘法完成，更准确地说是log2(b)次平方运算+至多log2(b)次乘法运算。你要说平方和乘法不都一样吗，No，平方要快些，a*a，只需要将a变换到频域，自乘，再变回来，2次FFT，而乘法是3次，快一半左右。

%]7q zU4R,f0除法怎么做？用乘法做，最简单的是用牛顿迭代，设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)要好。

t8CZLZq9j*Q0其实实现了乘法就差不多把整个高精做得差不多了，用你已经实现的再去实现另外其它的东西。你也可以先不去做除法，比如先做开方，或许有更好的方法回来再实现除法。

Gg(CGi!C^z0阶乘是几乎只依赖乘法的模块，做阶乘也挺有意思。阶乘从形式上看是连乘：

n!=1*2*...*n

zfl{z9}p0但是在着手进行连乘之间，希望做什么？使连乘的数变少一些，又想少做几次乘法，又想结果不出错？（鱼和熊掌？）51Testing软件测试网;r,x Q)A.}

&Oo _y'LgX0曾经一道ACM题是算阶乘末尾的0有多少，那个算法会启发你如何从阶乘中提取出素因子的幂，这样一来连乘的项数就会变少，但是遗憾的是如果n比较大素数也相当多，数论里一个经典的完美的证明，就是说这个世界素数是取之不尽数之不完。。。这里还牵涉到计算素数，因为阶乘增涨得比指数还快，数学分析里会让你算一个无穷小，分子是a^n，分母是n!，所以玩阶乘别玩上火，它发起脾气来上帝都受不了（不是只有上帝能追上指数吗），所以通常n较小，可n!很大，在2到n之间求素数，最方便的就是筛法，就是让所有正整数排个队，小的在前面，大的在后面，你从队伍里抓个最小的出来，很容易排头第一个就是，然后宣布他是素数，然后把所有是他的倍数的人踢出队伍，继续筛选，就得到全部的素数，当然当然，1是不算在内的，可以想像你就是1。

Hi6\ FZLe0阶乘变换成素因子的乘幂，再连乘，就是目前的1.0版，速度测试如下：

Factorial Function Test (1.0)
Y.~B7o8|p0_______________________
7fr.w,P-WzB `0n!      cost time (s)
$ckW(x}k}aKo{3@01k!     0.028163
F3VdO5Qz'S$C!hy010k!    0.856568
,FX:voqj020k!    2.215239
40k!    6.278020
80k!    16.57405851Testing软件测试网L'k z-F+Q)SH.@
Test Date\2007\04\11

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

:|X3a]:kwy{0Fibonacci Function Test51Testing软件测试网d3B{A#`%U E&X
_______________________
k-VZ%Is$G0f(n) cost time (s)
f(1k) 0.004677
b6p g\hb0f(10k) 0.058493
f(20k) 0.150679
f(40k) 0.286142
f(80k) 0.749612
T4|$ao r0

D2o+}E)HED?)k4?a#|0代码现在还在改进加测试（一砣砣的注释代码不敢轻易删），暂不提供。用于测试的代码是：

#include "TestTime.h"
)c'oMa&~wu#m|0TestTime tt;
${'V-^M$P(m"uv[ eo0int main(int argc, char* argv[])
5c2\Mm)dDbT%c$q0{51Testing软件测试网a:g6I&]2kw
//char *str=new char[200001];
+\S[B0c5L"bub"oM0 //unsigned long n;
0k'x3dAa0 hreal a(0);
{h&b8|0q0 double t1k,t10k,t20k,t40k,t80k;
8|+y8CN(Y^c } x*y0 //printf("n=");
9_0PC ^-s ZDfL+a[1R0 //scanf("%u",&n);
t1k=tt.currTime();
k:ovBr0 a.fibonacci(1000);
t1k=tt.currTime()-t1k;

t10k=tt.currTime();
~^cm!lJy0 a.fibonacci(10000);
t10k=tt.currTime()-t10k;

t20k=tt.currTime();51Testing软件测试网z7M'sdQ$T
a.fibonacci(20000);51Testing软件测试网/N!z]l JYf
t20k=tt.currTime()-t20k;51Testing软件测试网,C2QBu0i6b6n

wR!Xxp0 t40k=tt.currTime();51Testing软件测试网$Ky)C9_Qf0R
a.fibonacci(40000);
E8v,Y/N"_0 t40k=tt.currTime()-t40k;

#p,Gqk`0 t80k=tt.currTime();
L b"W*V;Z8o6~0 a.fibonacci(80000);51Testing软件测试网HJ%`Rt:U
t80k=tt.currTime()-t80k;

0l[}{ au0 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);
@G/YIM v A'w9n0 //a.display(str);
| v|F6C Ph Y0 //printf("%s\n",str);51Testing软件测试网D']6gX8i8{d
//printf("%u! finished.\ncost %lf(ms)\n",n,(t2-t1)*1000);
kY:z^cD%kEZ0 //delete[] str;
qwTY!p|0 //system("pause");51Testing软件测试网(Y nc YSA] Y
return 0;
b K#bc1P.aZ0}

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

!b)?F)q}w0#include <windows.h>

Qe/`2K/^8i(^x6J,d0class TestTime
{
lpB*GO+p0?nq0 double freq;
Gvu|"x0public:
$eWF0v``0z_k0 TestTime()
{
QVJ9Z}:Y N$t0  LARGE_INTEGER li_freq;51Testing软件测试网Y&`,s1P\i]
  if(!QueryPerformanceFrequency(&li_freq))51Testing软件测试网?'u#jo!xt
   freq=-1;
"q*_qcz0  else
HoW-M4f5`6N0  {
.G}aZ&i/?0   freq=li_freq.HighPart * 4294967296.0 + li_freq.LowPart;
3b]2`S4D%MW i'e0  }51Testing软件测试网/w(v#Q~5|&CG
}
~Ns"o9JM:]0 double currTime()
T \:R'v(["K0 {
.I\7yc7Hu+d0  LARGE_INTEGER nowCount;
  QueryPerformanceCounter(&nowCount);51Testing软件测试网1Lie5ui6jj
  double time=nowCount.HighPart * 4294967296.0 + nowCount.LowPart;
,S#I3{ S b~t(uGUI0  return time/freq;51Testing软件测试网,MW^l's\I
}
m+MwM"PX0};