FFT模板 - #include <Seter> - using namespace Orz;
FFT模板
坑了个爹的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; }
2011年10月22日 08:52
- -
2011年10月22日 09:28
被FFT虐爆了
2011年10月22日 09:51
只手动FT过的飘过……-_-;;
2011年10月22日 10:57
不会复数的路过……
2011年10月22日 11:26
- - 那个对你来说小case啦