FFT模板 - #include <Seter> - using namespace Orz;

FFT模板

Seter posted @ 2011年10月21日 23:07 in Practice with tags HighPrecision FFT , 3899 阅读

坑了个爹的FFT……原理只了解到卷积那一步……怎么转换成整系数而非sincos乱搞的系数那个完全不了解……权当背代码了……

其实代码灰常好写……根本不像我想像的那样超级复杂!代码的关键在于如何选P([tex]$P=C*2^{\log n}+1$[/tex])……

其实P是很好选的,对于一个n,用各种c测试下是不是素数就可以了……一般将P固定起来。我用的是[tex]$P=479*2^{21}+1$[/tex]。为什么用21位呢?再大一点就要unsigned了,不方便(尽管GAOA上建议使用[tex]$P=15*2^{27}+1$[/tex])。小一点又会过于小。

选了P以后要计算它的原根……这个就有点晕了,不过我可以告诉你,479<<21^1的原根是3,15<<27^1的原根是31……背出来就OK了。

由于我自己也不是完全懂……所以不细讲了……

蝴蝶操作:[tex]$x_{i+j}=x_{i+j}+G^{\frac{j(P-1)}{2^m}} x_{i+j+l}\,,\,x_{i+j+l}=x_{i+j}-G^{\frac{j(P-1)}{2^m}} x_{i+j+l}$[/tex]。原根的幂可以先预处理出来。

具体操作是先对AB进行FFT,再让B成为卷积之积,然后映射回A上做IDFT。由于我们换上了P,所以特别方便,只要把B倒着放回A里面,然后直接对A做FFT,除以n就是乘上[tex]$n^{P-2}\pmod{P}$[/tex]!神奇吧!

对于迭代时那个奇葩的04261537……最简单的方法是写个rev函数将i的二进制位倒过来……时间复杂度是O(lgn)的。GAOA上说有更好的方法,但没说,于是我YY了一个。注意到:

4=0+4

2=0+2,6=4+2

1=0+1,5=4+1,3=2+1,7=6+1

发现好像可以a[0]=x[0]后,先推1个,然后用前两个推下两个,再用前四个推下四个……就解决了。

总的要把握宁滥勿缺的原则……那啥……搞大点总不会错的嘛。

下面是代码……只有31行……不过估计没人看得懂 - -

#include <stdio.h>
#define maxn 1<<10
#define getstr(s,l) while((s[l]=getchar()-48)>=0&&s[l++]<=9);
#define P (479<<21^1)
int n,a[maxn<<1],b[maxn<<1],p[maxn],s,sa[maxn],sb[maxn];
long long g[21];
inline pw(long long x,unsigned y){
	long long ans=1;
	for(y<<=1;y>>=1;x=x*x%P)if(y&1)ans=ans*x%P;return ans;
}inline void mk(int *x,int *y){
	int i,j=1,l=1,of=n>>1;
	for(i=0;(i+=j)!=n;l<<=1,of>>=1)for(j=-1;++j!=l;)if((p[i+j]=p[i+j-l]-of)>=0)x[i+j]=y[p[i+j]];
}inline void FFT(int *x){
	int m,i,j,l=1,tk,tx;long long t;
	for(m=-1;++m!=s;l<<=1)for(i=-l<<1;(i+=l<<1)!=n;)for(j=-1,t=1;++j!=l;t=g[m]*t%P){
		x[i+j]=((tx=x[i+j])+(tk=t*x[i+j+l]%P))%P;
		if((x[i+j+l]=(tx-tk)%P)<0)x[i+j+l]+=P;
	}
}main(){
	int i,la=0,lb=0,t=479<<21;long long dft;
	getstr(sa,la);getstr(sb,lb);
	dft=pw(n=1<<(s=33-__builtin_clz((la>lb?la:lb)-1)),P-2);
	for(i=-1;++i!=s;g[i]=pw(3,t>>=1));
	a[0]=sa[p[0]=la-1];mk(a,sa);FFT(a);
	b[0]=sb[p[0]=lb-1];mk(b,sb);FFT(b);
	for(i=-1;++i!=n;b[i]=(long long)a[i]*b[i]%P);
	a[0]=b[0];p[0]=n;mk(a,b);FFT(a);
	for(a[la=i=0]=dft*a[0]%P;++i!=n;a[i-1]%=10)if(a[i]=dft*a[i]%P+a[i-1]/10)la=i;
	for(la++;la--;putchar(a[la]+48));putchar('\n');
	return 0;
}
By Seter
Avatar_small
λ 说:
2011年10月22日 09:51

只手动FT过的飘过……-_-;;

Avatar_small
Seter 说:
2011年10月22日 10:57

不会复数的路过……

Avatar_small
λ 说:
2011年10月22日 11:26

- - 那个对你来说小case啦


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter
Host by is-Programmer.com | Power by Chito 1.3.3 beta | © 2007 LinuxGem | Design by Matthew "Agent Spork" McGee