首页 > itarticle > lg1919 FFT NTT

lg1919 FFT NTT

admin 11月 13, 2020 0

链接:lg1919

简单FFT模板，将n位数分解为$a_0+a_1 * x+a_2 * x^2+...+a_{n-1} * x^{n-1}$,
$a_j=frac{1}{n} sum_{k=0}^{n-1} y_komega_n^{-kj}$
$$dp[i]=i^n-sum_{j=1}^{i-1}{C_i^j*dp[j]}$$

FFT


#define rep(i,x,y) for(register int i=x;i<=y;++i)
#define repd(i,x,y) for(register int i=x;i>=y;--i)
#define PI acos(-1)
using namespace std;
const int N=1e6+7;
struct {
	double x,y;
	cpx(){};
	cpx(double a,double b){x=a;y=b;}
	inline cpx operator * (cpx b){return cpx(x*b.x-y*b.y,b.x*y+b.y*x);}
	inline cpx operator *= (cpx b){*this=*this*b;}
	inline cpx operator + (cpx b){return cpx(x+b.x,y+b.y);}
	inline cpx operator - (cpx b){return cpx(x-b.x,y-b.y);}
}A[N],B[N];
char str[N];
int n,m,R[N],L,ans[N];
void fft(cpx *a,int f){
	rep(i,0,n)if(i<R[i])swap(a[i],a[R[i]]);
	for(register int i=1;i<n;i<<=1){
		cpx nw=cpx(cos(PI/i),f*sin(PI/i));
		for(register int j=0;j<n;j+=(i<<1)){
			cpx w=cpx(1,0);
			for(register int k=0;k<i;k++,w*=nw){
				cpx x=a[j+k],y=w*a[i+j+k];
				a[j+k]=x+y;a[i+j+k]=x-y;
			}
		}
	}
	if(f==-1)rep(i,0,n-1)a[i].x/=n;
}
int main(){
	scanf("%d",&n);
	scanf("%s",str);
	repd(i,n-1,0)A[n-1-i].x=str[i]-'0';
	scanf("%s",str);
	repd(i,n-1,0)B[n-1-i].x=str[i]-'0';
	m=n<<1;
	for(n=1;n<=m;n<<=1,++L);
	rep(i,0,n-1)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	fft(A,1);fft(B,1);
	rep(i,0,n-1)A[i]*=B[i];
	fft(A,-1);
	rep(i,0,n-1){
		ans[i]+=(int)(A[i].x+0.5);
		if(ans[i]>=10){
			ans[i+1]+=ans[i]/10;
			ans[i]%=10;
			n+=(i==n);
		}
	}
	while(!ans[n]&&n>=1)n--;
	m++;
	repd(i,n,0)printf("%d",ans[i]);
	return 0;
}

NTT


#define rep(i,x,y) for(register int i=x;i<=y;++i)
#define repd(i,x,y) for(register int i=x;i>=y;--i)
#define ll long long
#define mod 998244353
using namespace std;
const int N=1e6+7;
ll A[N],B[N],R[N];
int n,m,bit;
char str[N];
ll qpow(ll x,ll y){
	ll z=x,ans=1;
	for(;y;y>>=1,z=z*z%mod)if(y&1)ans=ans*z%mod;
	return ans%mod;
}
void ntt(ll *a,int f){
	rep(i,0,n-1)if(i<R[i])swap(a[i],a[R[i]]);
	for(int i=1;i<n;i<<=1){
		ll wn=qpow(3,(mod-1)/(i<<1));
		if(f==-1)wn=qpow(wn,mod-2);
		for(int j=0;j<n;j+=i<<1){
			ll w=1;
			for(int k=0;k<i;k++,w=w*wn%mod){
				ll x=a[j+k],y=w*a[i+j+k]%mod;
				a[j+k]=(x+y+mod)%mod;
				a[i+j+k]=(x-y+mod)%mod;
			}
		}
	}
	if(f==-1){
		ll inv=qpow(n,mod-2);
		rep(i,0,n-1)a[i]=a[i]*inv%mod;
	}
}
int main(){
	scanf("%d",&n);
	scanf("%s",str);
	repd(i,n-1,0)A[n-1-i]=str[i]-'0';
	scanf("%s",str);
	repd(i,n-1,0)B[n-1-i]=str[i]-'0';
	m=n<<1;
	for(n=1;n<=m;n<<=1,++bit);
	rep(i,0,n-1)R[i]=(R[i>>1]>>1)|((i&1)<<(bit-1));
	ntt(A,1);ntt(B,1);
	rep(i,0,n-1)A[i]=A[i]*B[i]%mod;
	ntt(A,-1);
	rep(i,0,n-1){
		if(A[i]>=10){
			A[i+1]+=A[i]/10;
			A[i]%=10;
		}
	}
	while(A[n-1]==0)--n;
	repd(i,n-1,0)printf("%d",A[i]);
	return 0;
}