【模板】拆系数FFT

#include<stdio.h>
#include<cmath>
char buf[1<<20],*p1,*p2;
#define GC (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<20,stdin),p1==p2)?0:*p1++)
inline int R(){
	register char t=GC;
	register int x=0;
	while(t<'!') t=GC;
	while(t>='!') x=x*10+t-48,t=GC;
	return x;
}
char pbuf[1<<20],*pp=pbuf;
inline void push(const char c){
	if(pp-pbuf==(1<<20))fwrite(pbuf,1,(1<<20),stdout),pp=pbuf;
	*pp++=c;
}
inline void write(int x){
	static int sta[35];
	register int top=0,tmp;
	do{tmp=x/10,sta[top++]=x-tmp*10,x=tmp;}while(x);
	while(top) push(sta[--top]+'0');
}
std::queue<int>ans;
const double pi=acos((long double)-1);
int inf;
struct comp{
	double a,b;
	comp(){}
	comp(double A,double B){
		a=A;b=B;
	}
	comp(double P){
		a=cos(P);b=sin(P);
	}
	inline comp operator+(const double &A)const{return comp(a+A,b);}
	inline comp operator+(const comp &A)const{return comp(a+A.a,b+A.b);}
	inline comp operator-(const double &A)const{return comp(a-A,b);}
	inline comp operator-(const comp &A)const{return comp(a-A.a,b-A.b);}
	inline comp operator*(const double &A)const{return comp(a*A,b*A);}
	inline comp operator*(const comp &A)const{return comp(a*A.a-b*A.b,a*A.b+b*A.a);}
}tmp1[1048576],tmp2[1048576],tmp3[1048576],tmp4[1048576],tmp5[1048576],tmp[1048576];
inline void Swap(comp &a,comp &b){comp tmp=a;a=b;b=tmp;}
int _0[1048576],_lset=-1,inv[1048576];
inline void mod(int &p){p-=p>=inf?inf:0;}
inline int ksm(long long a,int b){register int ans=1;while(b)(b&1)&&(ans=a*ans%inf),a=a*a%inf,b>>=1;return ans;}
inline void init(int n){if(_lset==n)return;_lset=n;_0[0]=0;for(register int i=1;i<n;i++)_0[i]=i&1?_0[i^1]|n>>1:_0[i>>1]>>1;}
inline void clr(int a[],int l,int r){for(register int i=l;i<r;i++)a[i]=0;}
inline void cpy(int a[],int b[],int n){for(register int i=0;i<n;i++)b[i]=a[i];}
inline int gett(int n){while(n!=(n&-n))n^=n&-n;return n<<1;}
inline int getint(double p){
	p+=0.5;
	return int(p-double(1)*inf*floor(p/double(inf)));
}
inline void dft(comp a[],int n,bool typ){
	init(n);register comp q;tmp[0]=comp(1,0);
	for(register int i=0;i<n;i++)if(i<_0[i])Swap(a[i],a[_0[i]]);
	for(register int i=1;i<n;i<<=1){
		comp w=comp(typ?-pi/i:pi/i);
		for(register int k=i-2;k>=0;k-=2)
			tmp[k^1]=(tmp[k]=tmp[k>>1])*w;
		for(register int j=0;j<n;j+=i<<1)
			for(register int k=0;k<i;k++)
				q=tmp[k]*a[i+j+k],a[i+j+k]=a[j+k]-q,a[j+k]=a[j+k]+q;
	}if(typ)for(register int i=0;i<n;i++)a[i].a/=double(n),a[i].b/=double(n);
}
inline void fft(int a[],int b[],int c[],int n){
	for(register int i=0;i<n;i++)tmp1[i]=comp(a[i]&32767,a[i]>>15),tmp2[i]=comp(b[i]&32767,b[i]>>15);
	dft(tmp1,n,0);dft(tmp2,n,0);
	for(register int i=0,j;i<n;i++)
		j=i?n-i:0,tmp3[i]=comp(tmp1[i].a+tmp1[j].a,tmp1[i].b-tmp1[j].b)*tmp2[i]*0.5,tmp4[i]=comp(tmp1[i].b+tmp1[j].b,tmp1[j].a-tmp1[i].a)*tmp2[i]*0.5;
	dft(tmp3,n,1);dft(tmp4,n,1);
	int val2=32768%inf,val3=1073741824%inf;
	for(register int i=0;i<n;i++){
		register int v1=getint(tmp3[i].a),v2=getint(tmp3[i].b+tmp4[i].a),v3=getint(tmp4[i].b);
		c[i]=(1ull*val2*v2+1ull*val3*v3+v1)%inf;
	}
}
inline void getinv(int a[],int b[],int tmp[],int n){
	for(register int i=0;i<n<<1;i++)b[i]=0;
	b[0]=ksm(a[0],inf-2);
	for(register int i=1;i<n;i<<=1){
		cpy(a,tmp,i<<1);
		fft(b,tmp,tmp,i<<2);
		mod(tmp[0]+=inf-2);
		for(register int j=0;j<i<<2;j++)mod(tmp[j]=inf-tmp[j]);
		fft(b,tmp,b,i<<2);
		clr(b,i<<1,i<<2);
	}
}