FFT

一直挺想学FFT。

1.n次单位根

http://wenku.baidu.com/link?url=x4Ot-WwMJOUCzs8soBnlVTQzlFJy0evAuVJYVG4Pe68ASMgws9LR1r7i02QZnx94JDSyoHNSb0s4dFObMyeTClbiNssF2dFqg5GZ0NVATJ_

2.FFT

http://blog.csdn.net/iamzky/article/details/22712347

然后。。。我发现。。。在网上找的资料我都看不懂。。。然后就看算法导论发现讲的非常详细。。。然后就看zky神犇的代码看看看看看看看看看看看看。。。。终于看懂了。。。据说还有NTT，CRT什么的。。。给我点智商行不行啊TAT。。。

以下是zky神犇的代码。

#include<cstdio>  
#include<cmath>  
#include<cstring>  
#include<iostream>  
#include<algorithm>  
using namespace std;  
const int maxn=1e6+10;  
struct cp{  
    double r,i;  
    cp(double _r=0,double _i=0):  
        r(_r),i(_i){}  
    cp operator+(cp x){return cp(r+x.r,i+x.i);}  
    cp operator-(cp x){return cp(r-x.r,i-x.i);}  
    cp operator*(cp x){return cp(r*x.r-i*x.i,r*x.i+i*x.r);}  
};  
cp a[maxn],b[maxn],A[maxn],x,y,c[maxn];  
char s1[maxn],s2[maxn];  
int sum[maxn],a1[maxn],a2[maxn],dig[maxn];  
int len1,len2,rev[maxn],N,L;  
void FFT(cp a[],int flag){  
    for(int i=0;i<N;i++)A[i]=a[rev[i]];  
    for(int i=0;i<N;i++)a[i]=A[i];  
    for(int i=2;i<=N;i<<=1){  
        cp wn(cos(2*M_PI/i),flag*sin(2*M_PI/i));  
        for(int k=0;k<N;k+=i){  
            cp w(1,0);  
            for(int j=k;j<k+i/2;j++){  
                x=a[j];  
                y=a[j+i/2]*w;  
                a[j]=x+y;  
                a[j+i/2]=x-y;  
                w=w*wn;  
            }  
        }  
    }  
    if(flag==-1)for(int i=0;i<N;i++)a[i].r/=N;  
}  
int main(){  
    scanf("%s%s",s1,s2);  
    len1=strlen(s1);  
    len2=strlen(s2);  
    for(N=1,L=0;N<max(len1,len2);N<<=1,L++);N<<=1;L++;  
    for(int i=0;i<N;i++){  
        int len=0;  
        for(int t=i;t;t>>=1)dig[len++]=t&1;  
        for(int j=0;j<L;j++)rev[i]=(rev[i]<<1)|dig[j];  
    }  
    for(int i=0;i<len1;i++)a1[len1-i-1]=s1[i]-'0';  
    for(int i=0;i<len2;i++)a2[len2-i-1]=s2[i]-'0';  
    for(int i=0;i<N;i++)a[i]=cp(a1[i]);  
    for(int i=0;i<N;i++)b[i]=cp(a2[i]);  
    FFT(a,1);FFT(b,1);  
    for(int i=0;i<N;i++)c[i]=a[i]*b[i];  
    FFT(c,-1);  
    for(int i=0;i<N;i++)sum[i]=c[i].r+0.5;  
    for(int i=0;i<N;i++){  
        sum[i+1]+=sum[i]/10;  
        sum[i]%=10;  
    }  
    int l=len1+len2-1;  
    while(sum[l]==0&&l>0)l--;  
    for(int i=l;i>=0;i--)  
    putchar(sum[i]+'0');  
    putchar('
');  
    return 0;  
}  

ps：http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform 讲的也挺详细了。。。后面的扩展还没有看懂TAT

ps：NTT：快速数论变换：http://blog.csdn.net/acdreamers/article/details/39026505 

ps：CRT：中国剩余定理: http://blog.csdn.net/acdreamers/article/details/8050018 

51nod1028：修改后的FFT

#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#include<cmath>
using namespace std;
#define rep(i,s,t) for(int i=s;i<=t;i++)
#define dwn(i,s,t) for(int i=s;i>=t;i--)
#define clr(x,c) memset(x,c,sizeof(x))
int read(){
	int x=0;char c=getchar();
	while(!isdigit(c)) c=getchar();
	while(isdigit(c)) x=x*10+c-'0',c=getchar();
	return x;
} 
const int nmax=4e5+5;
const int inf=0x7f7f7f7f;
struct c{
	double r,i;
	c(double r=0,double i=0):r(r),i(i){}
	c operator+(c x){
		return c(r+x.r,i+x.i);
	}
	c operator-(c x){
		return c(r-x.r,i-x.i);
	}
	c operator*(c x){
		return c(r*x.r-i*x.i,r*x.i+i*x.r);
	}
};
c a[nmax],b[nmax];int rev[nmax],ans[nmax],n;char sa[nmax],sb[nmax];
void fft(c *a,int f){
	rep(i,0,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=2;i<=n;i<<=1){
		c wn(cos(2*M_PI/i),f*sin(2*M_PI/i));
		for(int j=0;j<n;j+=i){
			c w(1,0);
			for(int k=0;k<i/2;++k,w=w*wn){
				c x=a[j+k],y=w*a[j+k+i/2];
				a[j+k]=x+y;a[j+k+i/2]=x-y;
			}
		}
	}
	if(f<0) rep(i,0,n-1) a[i].r/=n;
}
int main(){
	scanf("%s%s",sa,sb);
	int lena=strlen(sa)-1,lenb=strlen(sb)-1;n=max(lena,lenb);
	rep(i,0,lena) a[i]=c(sa[lena-i]-'0');
	rep(i,0,lenb) b[i]=c(sb[lenb-i]-'0');
	int m=n+n+2,L=0;for(n=1;n<=m;n<<=1,++L);//n是取不到的。那么1<<L就是取不到的。那么下一行应该是1<<L-1 
	rep(i,0,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
	fft(a,1);fft(b,1);
	rep(i,0,n-1) a[i]=a[i]*b[i];
	fft(a,-1);
	rep(i,0,m) ans[i]=a[i].r+0.1;
	rep(i,0,m) ans[i+1]+=ans[i]/10,ans[i]%=10;
	while(ans[m]==0&&m) --m;
	dwn(i,m,0) putchar(ans[i]+'0');putchar('
');
	return 0;
} 

bzoj2194：打一打模板

n=8 c[6]=a[7]b[1]+a[8]b[2]=a[7]b[7]+a[8]b[6] c[7]=a[8]b[1]...=a[8]b[7].所以将b逆向存一下就好了。速度进第一页啦嘿嘿嘿

#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#include<cmath>
using namespace std;
#define rep(i,s,t) for(int i=s;i<=t;i++)
#define dwn(i,s,t) for(int i=s;i>=t;i--)
#define clr(x,c) memset(x,c,sizeof(x))
int read(){
	int x=0;char c=getchar();
	while(!isdigit(c)) c=getchar();
	while(isdigit(c)) x=x*10+c-'0',c=getchar();
	return x;
} 
const int nmax=4e5+5;
const int inf=0x7f7f7f7f;
struct c{
	double r,i;
	c(double r=0,double i=0):r(r),i(i){};
	c operator+(c x){
		return c(r+x.r,i+x.i);
	}
	c operator-(c x){
		return c(r-x.r,i-x.i);
	}
	c operator*(c x){
		return c(r*x.r-i*x.i,r*x.i+i*x.r);
	}
};
c a[nmax],b[nmax];
int rev[nmax],ans[nmax],n;
void fft(c *a,int f){
	rep(i,0,n-1) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=2;i<=n;i<<=1){
		c wn(cos(2*M_PI/i),f*sin(2*M_PI/i));
		for(int j=0;j<n;j+=i){
			c w(1,0);
			for(int k=0;k<i/2;k++,w=w*wn){
				c x=a[j+k],y=w*a[j+k+i/2];
				a[j+k]=x+y;a[j+k+i/2]=x-y;
			}
		}
	}
	if(f<0) rep(i,0,n-1) a[i].r/=n;
}
int main(){
	n=read();--n;
	rep(i,0,n) a[i]=c(read()),b[n-i]=c(read());
	int m=n+n,L=0;for(n=1;n<=m;n<<=1,L++);
	rep(i,0,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
	fft(a,1);fft(b,1);
	rep(i,0,n-1) a[i]=a[i]*b[i];
	fft(a,-1);
	rep(i,m/2,m) printf("%d
",(int)(a[i].r+0.1));
	return 0;
}
/*
5
3 1
2 4
1 1
2 4
1 4
*/

bzoj3527：力

话说还是不懂卷积是什么东西啊TAT。以下是我网上找到的神犇详细题解 by晴歌 http://www.cnblogs.com/Sunnie69/p/5574222.html

其实把m设成n+n+10一般就不会出现边界问题还有<=n不要是<n就可以了。话说这种写法优美多了

#include<cstdio>
#include<cstring>
#include<cctype>
#include<algorithm>
#include<cmath>
using namespace std;
#define rep(i,s,t) for(int i=s;i<=t;i++)
#define dwn(i,s,t) for(int i=s;i>=t;i--)
#define clr(x,c) memset(x,c,sizeof(x))
int read(){
	int x=0;char c=getchar();
	while(!isdigit(c)) c=getchar();
	while(isdigit(c)) x=x*10+c-'0',c=getchar();
	return x;
}
const int nmax=3e5+5;//三倍空间就可以了 
const int inf=0x7f7f7f7f;
struct c{
	double r,i;
	c(double r=0,double i=0):r(r),i(i){}
	c operator+(c x){
		return c(r+x.r,i+x.i);
	}
	c operator-(c x){
		return c(r-x.r,i-x.i);
	}
	c operator*(c x){
		return c(r*x.r-i*x.i,r*x.i+i*x.r);
	}
};
c a[nmax],b[nmax];int rev[nmax],n;
double f[nmax],g[nmax],t[nmax],ans[nmax];
void dft(c *a,int f){
	rep(i,0,n) if(i<rev[i]) swap(a[i],a[rev[i]]);
	for(int i=2;i<=n;i<<=1){
		c wn=c(cos(2*M_PI/i),f*sin(2*M_PI/i)); //因为这里忘了f*一直调TAT 
		for(int j=0;j<n;j+=i){
			c w=c(1,0);
			for(int k=0;k<i/2;++k,w=w*wn){
				c x=a[j+k];c y=w*a[j+k+i/2];
				a[j+k]=x+y;a[j+k+i/2]=x-y;
			}
		}
	}
	if(f<0) rep(i,0,n) a[i].r/=n;
}
void fft(double *ta,double *tb,c *a,c *b){
	rep(i,0,n) a[i]=c(ta[i]),b[i]=c(tb[i]);
	dft(a,1);dft(b,1);
	rep(i,0,n) a[i]=a[i]*b[i];
	dft(a,-1);
}
int main(){
	n=read();
	rep(i,1,n) scanf("%lf",&f[i]),g[i]=1.0/i/i;
	int m=n+n+10,L=0,tm=n;for(n=1;n<=m;n<<=1,L++);
	rep(i,0,n) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
	fft(f,g,a,b);
	rep(i,0,tm) ans[i]=a[i].r,t[i]=f[tm-i];
	fft(t,g,a,b);
	rep(i,0,tm) ans[i]-=a[tm-i].r;
	rep(i,1,tm) printf("%.3lf
",ans[i]);
	return 0;
}
/*
5  
4006373.885184 
15375036.435759 
1717456.469144 
8514941.004912 
1410681.345880 
 */