[大坑]FFT学习

[大坑]FFT学习

Macros

#define fon(i,s)    for(int i=0;i<s; ++i)
#define fone(i,s)   for(int i=0;i<=s;++i)
#define fox(i,f,t)  for(int i=f;i<t; ++i)
#define foxe(i,f,t) for(int i=f;i<=t;++i)
#define don(i,s)    for(int i=s;i;   --i)
#define done(i,s)   for(int i=s;~i;  --i)
#define dox(i,f,t)  for(int i=f;i>t; --i)
#define doxe(i,f,t) for(int i=f;i>=t;--i)
#define ifm(a,b) if((a)<(b))
#define _swp(a,b) std::swap(a,b)
#define lp while(1)
#define qlp break;
#define nlp continue;
#define maxp 30
#define odd(x) (x&1)
#define even(x) !(x&1)
#define _cl(x,f,t) fox(_CLEAR,f,t) x[_CLEAR]=0
template<class T> inline void _st(T* f,T* t,T p){
	for(T* x=f;x<t;++x) *x=p;
}

Bit Reverse

inline void _BR(int* a,int r){
	for(int i=0,j=1;i<r;++i,j<<=1){
		for(int k=0,kx=j;k<j;++k,++kx){
			a[k]=a[k]<<1;
			a[kx]=a[k]|1;
		}
	}
}
inline void _BR_iter(int* a,int r){
	int u=r;
	fon(i,r){
		a[i]=a[i]<<1;
		a[u++]=a[i]|1;
	}
}
inline void _BR_diter(int* a,int r){
	fon(i,r) a[i]>>=1;
}

Fast power mod

wjz大爷说他的fpm只要一行吓cry.
经典沙茶zbt写法.

inline int fpm(int a,int b,int p){
	int q=1;
	while(b){
		if(b&1) q=((long long)q*a)%p;
		a=((long long)a*a)%p;
		b>>=1;
	}
	return q;
}

NTT

感觉FFT和IFFT分开来写会好一些→ →

struct _NTT_base{
	int mod,w1,wm;
	int p[maxp],pi[maxp],d;
	inline int inv(int p){
		return fpm(p,mod-2,mod);
	}
	inline void init(int m,int w){
		mod=m,p[0]=w1=w;
		int u=m-1,u2=m-1;
		d=0;
		while(even(u2)) u2>>=1;
		p[0]=fpm(p[0],u2,m);
		pi[0]=inv(p[0]);
		while(even(u)){
			++d;
			p[d]=((long long)p[d-1]*p[d-1])%m,pi[d]=((long long)pi[d-1]*pi[d-1])%m;
			u>>=1;
		}
	}
	inline void FFT(int* a,int* bitrev,int l){
		fon(i,l) ifm(i,bitrev[i]) _swp(a[i],a[bitrev[i]]);
		for(int i=2,h=1,xn=d-1;i<=l;i<<=1,h<<=1,--xn){
			int u=p[xn];
			for(int j=0;j<l;j+=i){
				int w=1;
				fox(k,j,j+h){
					int A=a[k],B=(long long)a[k+h]*w%mod;
					a[k]=(A+B)%mod,a[k+h]=(A-B+mod)%mod;
					w=(long long)w*u%mod;
				}
			}
		}
	}
	inline void IFFT(int* a,int* bitrev,int l){
		fon(i,l) ifm(i,bitrev[i]) _swp(a[i],a[bitrev[i]]);
		int invA=1,invB=(mod+1)>>1,invC=0;
		for(int i=2,h=1,xn=d-1;i<=l;i<<=1,h<<=1,--xn){
			int u=pi[xn];
			invA=(long long)invB*invA%mod;
			for(int j=0;j<l;j+=i){
				int w=1;
				fox(k,j,j+h){
					int A=a[k],B=(long long)a[k+h]*w%mod;
					a[k]=(A+B)%mod,a[k+h]=(A-B+mod)%mod;
					w=(long long)w*u%mod;
				}
			}
		}
		fon(i,l) a[i]=(long long)a[i]*invA%mod;
	}
	inline void FFT(int* a,int* b,int* bitrev,int l){
		fon(i,l) ifm(i,bitrev[i]) _swp(a[i],a[bitrev[i]]),_swp(b[i],b[bitrev[i]]);
		for(int i=2,h=1,xn=d-1;i<=l;i<<=1,h<<=1,--xn){
			int u=p[xn];
			for(int j=0;j<l;j+=i){
				int w=1;
				fox(k,j,j+h){
					int A=a[k],C=b[k],B=(long long)a[k+h]*w%mod,D=(long long)b[k+h]*w%mod;
					a[k]=(A+B)%mod,a[k+h]=(A-B+mod)%mod,b[k]=(C+D)%mod,b[k+h]=(C-D+mod)%mod;
					w=(long long)w*u%mod;
				}
			}
		}
	}
	inline void IFFT(int* a,int* b,int* bitrev,int l){
		fon(i,l) ifm(i,bitrev[i]) _swp(a[i],a[bitrev[i]]),_swp(b[i],b[bitrev[i]]);
		int invA=1,invB=(mod+1)>>1;
		for(int i=2,h=1,xn=d-1;i<=l;i<<=1,h<<=1,--xn){
			int u=pi[xn];
			invA=(long long)invA*invB%mod;
			for(int j=0;j<l;j+=i){
				int w=1;
				fox(k,j,j+h){
					int A=a[k],C=b[k],B=(long long)a[k+h]*w%mod,D=(long long)b[k+h]*w%mod;
					a[k]=(A+B)%mod,a[k+h]=(A-B+mod)%mod,b[k]=(C+D)%mod,b[k+h]=(C-D+mod)%mod;
					w=(long long)w*u%mod;
				}
			}
		}
		fon(i,l) a[i]=(long long)a[i]*invA%mod,b[i]=(long long)b[i]*invA%mod;
	}
};

这个(K^{-1}mod P)求法比较诡异...先求出(2^{-1}mod P)就是(frac{P+1}{2})(这个非常显然> <,P得是(2^kcdot c+1)所以是奇数),然后倍增,由于(K=2^u)...为了更好地运用循环资源> >...

坑点笔记

• in fpm(): + b>>=1;
• in _NTT_base::init(): int d error -> d
• in _NTT_base::IFFT(): calc invA method + invA*=invB - invA=invB,invB=invB*invB