2017-05-05

【算法笔记】多项式快速算法——FFT/NTT的应用

本文为博主原创，未经许可不得转载

优化高精度乘法

朴素的高精度乘法的复杂度是$O(n^2)$的

应用FFT可以做到$O(n\log n)$

我们可以把数写成多项式的形式,例如：

$1523=3\times 10^0 + 2 \times 10^1+5\times 10^2+1\times 10^3$

这样两数相乘实际上就是多项式相乘，上FFT就行了

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<ctime>
#include<cmath>
#include<algorithm>
typedef long long ll;
#define FILE "read"
#define MAXN 300010
//#define up(i,j,n) for(int i=j;i<=n;i++)
//#define dn(i,j,n) for(int i=j;i>=n;i--)
const double pi=acos(-1);
struct complex{
	double r,v;
	complex(double a=0,double b=0):r(a),v(b){}
	inline complex operator+(const complex& b){return complex(r+b.r,v+b.v);}
	inline complex operator-(const complex& b){return complex(r-b.r,v-b.v);}
	inline complex operator*(const complex& b){return complex(r*b.r-v*b.v,v*b.r+r*b.v);}
}w[MAXN],a[MAXN],b[MAXN];
char buf[1<<15],*fs,*ft;
inline char getc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
inline int read(){
	int x=0,f=1;  char ch=getc();
	while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getc();}
	while(isdigit(ch))  {x=x*10+ch-'0';  ch=getc();}
	return x*f;
}
inline void swap(complex& a,complex& b){complex t(a);a=b;b=t;}
char ch1[MAXN],ch2[MAXN];
int n,m,H,L,R[MAXN],ans[MAXN];
void FFT(complex* a,int f){
	up(i,0,L-1)  if(i<R[i])  swap(a[i],a[R[i]]);
	for(int len=2;len<=L;len<<=1){
		int l=len>>1;  //complex wn(cos(pi/l),f*sin(pi/l));
		up(i,1,l-1)  w[i]=w[i-1]*wn;
		for(int st=0;st<L;st+=len) up(k,0,l-1){
			complex x=a[st+k],y=w[k]*a[st+k+len/2];
			a[st+k]=x+y;  a[st+k+len/2]=x-y;
		}
	}
}
int main(){
	scanf("%s%s",ch1+1,ch2+1); w[0].r=1;
	n=strlen(ch1+1);  m=strlen(ch2+1);
	for(int len=0,i=n;i;i--,len++)  a[len]=ch1[i]-'0';
	for(int len=0,i=m;i;i--,len++)  b[len]=ch2[i]-'0';
	for(L=1,H=0;L<n+m;L<<=1,H++);
	up(i,0,L-1)  R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
	FFT(a,1);  FFT(b,1);
	up(i,0,L-1)  a[i]=a[i]*b[i];
	FFT(a,-1);
	int len=300000;
	up(i,0,n+m-2)  ans[i]=ll(a[i].r/L+0.5);
	up(i,0,len-1)  ans[i+1]+=ans[i]/10,ans[i]%=10;
	while(ans[len-1]==0)  len--;
	dn(i,len-1,0)  printf("%d",ans[i]);
	return 0;
}

计算卷积

形如$C_k=\sum_{j=0}^{k} A_j B_{k-j}$的式子叫作卷积

利用FFT可以快速求出卷积

【bzoj2194】传送门

我们把数组来回拧一拧，就成了卷积形式，然后上FFT

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 265000
const double pi=acos(-1);
char buf[1<<15],*fs,*ft;
int n,R[MAXN],L,H;
struct complex{
    double r,v;
    complex(double a=0,double b=0):r(a),v(b){}
    inline complex operator+(complex &b){return complex(r+b.r,v+b.v);}
    inline complex operator-(complex &b){return complex(r-b.r,v-b.v);}
    inline complex operator*(complex &b){return complex(r*b.r-v*b.v,r*b.v+v*b.r);}
}a[MAXN],b[MAXN],w[MAXN];
inline void swap(complex &a,complex &b){complex t(a);a=b;b=t;}
inline char getc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
inline int read(){
    int x=0,f=1;  char ch=getc();
    while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getc();}
    while(isdigit(ch))  {x=x*10+ch-'0';  ch=getc();}
    return x*f;
}
void FFT(complex *a,int f){
    for(int i=0;i<L;++i)  if(i<R[i])  swap(a[i],a[R[i]]);
    for(int len=2;len<=L;len<<=1){
        int l=len>>1;  //complex wn(cos(pi/l),f*sin(pi/l));
        for(int i=1;i<l;++i)  w[i]=w[i-1]*wn;
        for(int st=0;st<L;st+=len) for(int k=0;k<l;++k){
            complex x=a[st+k],y=w[k]*a[st+k+l];
            a[st+k]=x+y;  a[st+k+l]=x-y;
        }
    }
}
void solve(){
    n=read();  w[0].r=1;
    for(int i=0;i<n;++i) a[i]=read(),b[n-i-1]=read();
    for(L=1,H=0;L<n+n;L<<=1,++H);
    for(int i=0;i<L;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
    FFT(a,1);  FFT(b,1);
    for(int i=0;i<L;++i) a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<n;++i) printf("%d\n",(int)(a[i+n-1].r/L+0.5));
}
int main(){
    solve();
    return 0;
}

统计和为某值的数的对数

例一

【HDU4609】3-idiots

AC通道

题目大意：给出n个正整数，求从中任选3个数(下标不重复)，求以这三个数为长度的边能构成三角形的概率。

分析：令a[i]表示数i有多少个，对a求一遍卷积就可以得到任选两个数，和为x的方案数，这里注意要去一下重

由于三角形两边之和大于第三边，我们就能先用前缀和算出组不成三角形的方案，然后减一下就好了

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 300100
typedef long long ll;
const double pi=acos(-1);
struct complex{
	double r,v;
	complex(double a=0,double b=0):r(a),v(b){}
	inline complex operator +(const complex &b){return complex(r+b.r,v+b.v);}
	inline complex operator -(const complex &b){return complex(r-b.r,v-b.v);}
	inline complex operator *(const complex &b){return complex(r*b.r-v*b.v,r*b.v+v*b.r);}
}a[MAXN],w[MAXN];
int n,v[MAXN],R[MAXN];
ll b[MAXN],sum[MAXN];
inline int read(){
	int x=0,f=1;  char ch=getchar();
	while(ch>'9'||ch<'0')  {if(ch=='-')  f=-1;  ch=getchar();}
	while(ch>='0'&&ch<='9')  {x=x*10+ch-'0';  ch=getchar();}
	return x*f;
}
inline void swap(complex &a,complex &b){complex t(a);a=b;b=t;}
void FFT(complex *a,int f,int L){
	for(int i=0;i<L;++i) if(i<R[i]) swap(a[i],a[R[i]]);
	for(int len=2;len<=L;len<<=1){
		int l=(len>>1); //complex wn(cos(pi/l),f*sin(pi/l));
		for(int i=1;i<l;++i) w[i]=w[i-1]*wn;
		for(int st=0;st<L;st+=len)for(int k=0;k<l;++k){
			complex x=a[st+k],y=w[k]*a[st+k+l];
			a[st+k]=x+y; a[st+k+l]=x-y;
		}
	}
	if(f==-1) for(int i=0;i<L;++i) a[i].r/=L;
}
void solve(){
	n=read(); int L=1,H=0; w[0]=1; ll ans=0;
	memset(a,0,sizeof(a));
	memset(v,0,sizeof(v));
	for(int i=1;i<=n;++i)v[i]=read();
	for(int i=1;i<=n;++i)a[v[i]].r++;
	while(L<=200000)L<<=1,++H;
	for(int i=0;i<L;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
	FFT(a,1,L);
	for(int i=0;i<L;++i)a[i]=a[i]*a[i];
	FFT(a,-1,L);
	for(int i=0;i<L;++i)b[i]=(ll)(a[i].r+0.5);
	for(int i=1;i<=n;++i)b[v[i]*2]--;
	for(int i=0;i<L;++i)b[i]>>=1;
	for(int i=1;i<L;++i)sum[i]=sum[i-1]+b[i];
	for(int i=1;i<=n;++i)ans+=sum[v[i]];
	ll temp=(ll)n*(n-1)*(n-2)/6;
	printf("%.7lf\n",(double)(temp-ans)/temp);
}
int main(){
	//freopen(FILE".in","r",stdin);
	//freopen(FILE".out","w",stdout);
	for(int T=read();T;--T) solve();
	return 0;
}

例二

【SPOJ-TSUM】Triple Sums

题目大意：给出n(n<=40000)个整数（绝对值<=20000)，对任意三个下标不同的数求和，问和为-60000到60000的数对各有多少个。

分析：这次变成了三个，唯一的不同点在于去重，用容斥原理搞一搞就行了

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 265000
const double pi=acos(-1);
int n,L,H,R[MAXN],v[MAXN],c[MAXN],cnt1[MAXN],cnt2[MAXN],ans[MAXN];
char buf[1<<15],*fs,*ft;
struct complex{
	double r,v;
	complex(double a=0,double b=0):r(a),v(b){}
	inline complex operator+(const complex &b){return complex(r+b.r,v+b.v);}
	inline complex operator-(const complex &b){return complex(r-b.r,v-b.v);}
	inline complex operator*(const complex &b){return complex(r*b.r-v*b.v,r*b.v+v*b.r);}
}a[MAXN],b[MAXN],w[MAXN];
inline void swap(complex &a,complex &b){complex t(a);a=b;b=t;}
inline char getc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
inline int read(){
	int x=0,f=1;  char ch=getc();
	while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getc();}
	while(isdigit(ch))  {x=x*10+ch-'0';  ch=getc();}
	return x*f;
}
void FFT(complex *a,int f){
	for(int i=0;i<L;++i)  if(i<R[i]) swap(a[i],a[R[i]]);
	for(int len=2;len<=L;len<<=1){
		int l=len>>1;  //complex wn(cos(pi/l),f*sin(pi/l));
		for(int i=1;i<l;++i)  w[i]=w[i-1]*wn;
		for(int st=0;st<L;st+=len) for(int k=0;k<l;++k){
			complex x=a[st+k],y=w[k]*a[st+k+l];
			a[st+k]=x+y;  a[st+k+l]=x-y;
		}
	}
}
void solve(){
	n=read(); w[0].r=1;
	for(int i=1;i<=n;++i) v[i]=read()+20000;
	for(int i=1;i<=n;++i) a[v[i]].r++,b[v[i]*2].r++,c[v[i]*3]++;
	for(L=1,H=0;L<120000;L<<=1,++H);
	for(int i=0;i<L;++i)  R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
	FFT(a,1);  FFT(b,1);
	for(int i=0;i<L;++i)  b[i]=b[i]*a[i],a[i]=a[i]*a[i]*a[i];
	FFT(a,-1);  FFT(b,-1);
	for(int i=0;i<L;++i)  cnt1[i]=(int)(a[i].r/L+0.5),cnt2[i]=(int)(b[i].r/L+0.5);
	for(int i=0;i<L;++i)  ans[i]=(cnt1[i]-cnt2[i]*3+c[i]*2)/6;
	for(int i=0;i<L;++i)  if(ans[i])  printf("%d : %d\n",i-60000,ans[i]);
}
int main(){
	solve();
	return 0;
}

例三

【bzoj3771】Triple

这道题就是把例一和例二加起来

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 265000
const double pi=acos(-1);
int n,L,H,R[MAXN],v[MAXN],c1[MAXN],c2[MAXN],c3[MAXN],ans[MAXN];
char buf[1<<15],*fs,*ft;
struct complex{
    double r,v;
    complex(double a=0,double b=0):r(a),v(b){}
    inline complex operator+(const complex &b){return complex(r+b.r,v+b.v);}
    inline complex operator-(const complex &b){return complex(r-b.r,v-b.v);}
    inline complex operator*(const complex &b){return complex(r*b.r-v*b.v,r*b.v+v*b.r);}
}w[MAXN],a[MAXN],b[MAXN],c[MAXN];
inline void swap(complex &a,complex &b){complex t(a);a=b;b=t;}
inline char getc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
inline int read(){
    int x=0,f=1;  char ch=getc();
    while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getc();}
    while(isdigit(ch))  {x=x*10+ch-'0';  ch=getc();}
    return x*f;
}
void FFT(complex *a,int f){
    for(int i=0;i<L;++i)  if(i<R[i])  swap(a[i],a[R[i]]);
    for(int len=2;len<=L;len<<=1){
        int l=len>>1;  //complex wn(cos(pi/l),f*sin(pi/l));
        for(int i=1;i<l;++i)  w[i]=w[i-1]*wn;
        for(int st=0;st<L;st+=len) for(int k=0;k<l;++k){
            complex x=a[st+k],y=w[k]*a[st+k+l];
            a[st+k]=x+y;  a[st+k+l]=x-y;
        }
    }
}
void solve(){
    n=read();  w[0].r=1;
    for(int i=1;i<=n;++i)  v[i]=read();
    for(int i=1;i<=n;++i)  c1[v[i]]++,c2[v[i]*2]++,c3[v[i]*3]++;
    for(L=1,H=0;L<120000;L<<=1,++H);
    for(int i=0;i<L;++i)  a[i].r=c1[i],b[i].r=c2[i];
    for(int i=0;i<L;++i)  R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
    FFT(a,1);  FFT(b,1);
    for(int i=0;i<L;++i)  b[i]=b[i]*a[i],c[i]=a[i]*a[i],a[i]=a[i]*a[i]*a[i];
    FFT(a,-1);  FFT(b,-1);  FFT(c,-1);
    for(int i=0;i<L;++i)  ans[i]=((int)(a[i].r/L+0.5)-(int)(b[i].r/L+0.5)*3+c3[i]*2)/6+c1[i]+((int)(c[i].r/L+0.5)-c2[i])/2;
    for(int i=0;i<L;++i)  if(ans[i])  printf("%d %d\n",i,ans[i]);
}
int main(){
    solve();
    return 0;
}

例四

【uva12633】Super Rooks on Chessboard

题目大意：一个n*m的地图上有一些格子里有骑士，每个骑士可以攻击到它所在行，所在列以及所在的从左上到右下的对角线的所有格子，问有多少格子没被任何一个骑士攻击到。

数据范围：n,m<=50000,骑士数<=50000。

分析：如果以左下角为原点，会发现一条对角线上的(x+y)坐标是相同的，且任意两条对角线的(x+y)不同。由此我们就可以联想到卷积了。

设两个数组a,b,a[i]表示第i行是否被占领，b[i]表示第i列是否被占领，是的话值为0，否则值为1。

两个数组求个卷积，把没被攻击的第(x+y）项系数统计出来就行了。

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 265000
typedef long long ll;
const double pi=acos(-1);
char buf[1<<15],*fs,*ft;
ll T,ans,RR,CC,n,L,H,R[MAXN],check[MAXN];
struct complex{
	double r,v;
	complex(double a=0,double b=0):r(a),v(b){}
	inline complex operator+(const complex &b){return complex(r+b.r,v+b.v);}
	inline complex operator-(const complex &b){return complex(r-b.r,v-b.v);}
	inline complex operator*(const complex &b){return complex(r*b.r-v*b.v,r*b.v+v*b.r);}
}w[MAXN],a[MAXN],b[MAXN];
inline void swap(complex &a,complex &b){complex t(a);a=b;b=t;}
inline char getc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
inline ll read(){
	ll x=0,f=1;  char ch=getc();
	while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getc();}
	while(isdigit(ch))  {x=x*10+ch-'0';  ch=getc();}
	return x*f;
}
void FFT(complex *a,ll f){
	for(ll i=0;i<L;++i) if(i<R[i]) swap(a[i],a[R[i]]);
	for(ll len=2;len<=L;len<<=1){
		ll l=len>>1; //complex wn(cos(pi/l),f*sin(pi/l));
		for(ll i=1;i<l;++i)  w[i]=w[i-1]*wn;
		for(ll st=0;st<L;st+=len) for(ll k=0;k<l;++k){
			complex x=a[st+k],y=w[k]*a[st+k+l];
			a[st+k]=x+y;  a[st+k+l]=x-y;
		}
	}
}
void solve(){
	RR=read(); CC=read(); n=read(); w[0].r=1; ans=0;
	memset(a,0,sizeof(a));  memset(b,0,sizeof(b));
	for(ll i=0;i<RR;++i)  a[i].r=1;
	for(ll i=0;i<CC;++i)  b[i].r=1;
	memset(check,0,sizeof(check));
	for(ll i=1;i<=n;++i){
		ll x=RR-read(),y=read()-1;
		a[x].r=0;  b[y].r=0;  check[x+y]=1;
	}
	for(L=1,H=0;L<RR+CC;L<<=1,++H);
	for(ll i=0;i<L;++i)  R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
	FFT(a,1);  FFT(b,1);
	for(ll i=0;i<L;++i)  a[i]=a[i]*b[i];
	FFT(a,-1);
	for(ll i=0;i<RR+CC-1;++i) if(!check[i]) ans+=(ll)(a[i].r/L+0.5);
	printf("%lld\n",ans);
}
int main(){
	T=read();
	for(ll i=1;i<=T;++i)
		printf("Case %lld: ",i),solve();
	return 0;
}

解决01串匹配问题

例一

【bzoj4503】两个串

分析：定义两个长度相同的子串距离为$\sum (a_i-b_i)^2 b_i$

通配符置为0，这样就可以保证距离为0的两个串可以匹配，然后把s2反转一下求卷积就可以了。

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 300100
const double pi=acos(-1);
int n,m,temp,R[MAXN],f[MAXN],ans[MAXN];
char a[MAXN],b[MAXN];
struct complex{
    double r,v;
    complex(double a=0,double b=0):r(a),v(b){}
    inline complex operator +(const complex &b){return complex(r+b.r,v+b.v);}
    inline complex operator -(const complex &b){return complex(r-b.r,v-b.v);}
    inline complex operator *(const complex &b){return complex(r*b.r-v*b.v,r*b.v+v*b.r);}
}A[MAXN],B[MAXN],C[MAXN],D[MAXN],w[MAXN];
inline int read(){
    int x=0,f=1;  char ch=getchar();
    while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getchar();}
    while(isdigit(ch))  {x=x*10+ch-'0';  ch=getchar();}
    return x*f;
}
inline void swap(complex &a,complex &b){complex t(a);a=b;b=t;}
void FFT(complex *a,int f,int L){
    for(int i=0;i<L;++i) if(i<R[i]) swap(a[i],a[R[i]]);
    for(int len=2;len<=L;len<<=1){
        int l=(len>>1);  //complex wn(cos(pi/l),f*sin(pi/l));
        for(int i=1;i<l;++i)  w[i]=w[i-1]*wn;
        for(int st=0;st<L;st+=len)for(int k=0;k<l;++k){
            complex x=a[st+k],y=w[k]*a[st+k+l];
            a[st+k]=x+y;  a[st+k+l]=x-y;
        }
    }
    if(f==-1) for(int i=0;i<L;++i) a[i].r/=L;
}
int main(){
    //freopen(FILE".in","r",stdin);
    //freopen(FILE".out","w",stdout);
    scanf("%s%s",a,b);
    n=strlen(a);  m=strlen(b);
    for(int i=0;i<n;++i) A[i]=a[i]-'a'+1;
    for(int i=0;i<m;++i) B[i]=(b[i]=='?'?0:b[i]-'a'+1);
    for(int i=0;i<(m>>1);++i) swap(B[i],B[m-i-1]);
    for(int i=0;i<n;++i) C[i]=A[i]*A[i];
    for(int i=0;i<m;++i) D[i]=B[i]*B[i];
    for(int i=0;i<m;++i) temp+=B[i].r*B[i].r*B[i].r;
    for(int i=0;i<n;++i) A[i]=A[i]*2;
    int L=1,H=0; w[0]=1;
    while(L<=n+m) L<<=1,++H;
    for(int i=0;i<L;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
    FFT(A,1,L);  FFT(B,1,L);  FFT(C,1,L);  FFT(D,1,L);
    for(int i=0;i<L;++i) A[i]=C[i]*B[i]-A[i]*D[i];
    FFT(A,-1,L);
    for(int i=0;i<L;++i) f[i]=(int)(A[i].r+temp+0.5);
    for(int i=0;i<=n-m;++i) if(f[m-1+i]==0) ans[++ans[0]]=i;
    printf("%d\n",ans[0]);
    for(int i=1;i<=ans[0];++i) printf("%d\n",ans[i]);
    return 0;
}

优化dp

例一

【ZOJ3874】Permutation Graph

题目大意：给出一个1-n的排列，对于任意一个逆序对，在这两个位置连一条无向边，给你最终每个连通块的元素，求有多少种排列方案。

分析：很容易发现，每个连通块的元素一定是连续的，并且连通块之间顺序是固定的。那么我们可以预处理出长度为i的连通块的方案数，求解复杂度就可以和输入同阶了。

令dp[i]表示长度为i的连通块方案数。

Dp[i]=i!−∑𝑘!∗𝑑𝑝[i-k] (k=1…i),dp[0]=0.

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 300100
#define INF 1e9
#define cmin(a,b) a=min(a,b)
#define cmax(a,b) a=max(a,b)
using namespace std;
const int mod=786433;
const int G=10;
int n,m,fac[MAXN],dp[MAXN],a[MAXN],b[MAXN],R[MAXN],w[MAXN],wn[2][20];
inline int read(){
	int x=0,f=1;  char ch=getchar();
	while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getchar();}
	while(isdigit(ch))  {x=x*10+ch-'0';  ch=getchar();}
	return x*f;
}
inline void swap(int &a,int &b){a=a^b,b=a^b,a=a^b;}
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline int sub(int a,int b){return (a-=b)<0?a+mod:a;}
inline int mul(int a,int b){return 1LL*a*b%mod;}
inline int fast(int a,int b){int ret(1);for(;b;b>>=1,a=mul(a,a))if(b&1)ret=mul(ret,a);return ret;}
void NTT(int *a,int H,int f=0){
	int L=(1<<H);
	for(int i=0;i<L;++i) if(i<R[i]) swap(a[i],a[R[i]]);
	for(int e=1;e<=H;++e){
		int len=(1<<e),l=(len>>1);
		for(int i=1;i<l;++i) w[i]=mul(w[i-1],wn[f][e]);
		for(int st=0;st<L;st+=len)for(int k=0;k<l;++k){
			int x=a[st+k],y=mul(w[k],a[st+k+l]);
			a[st+k]=add(x,y); a[st+k+l]=sub(x,y);
		}
	}
	if(f) for(int i=0;i<L;++i) a[i]=mul(a[i],fast(L,mod-2));
}
void CDQ(int l,int r){
	if(l==r)  return;
	int mid=(l+r)>>1;
	CDQ(l,mid); 
	int L=1,H=0; while(L<=(r-l+1)) L<<=1,++H;
	for(int i=1;i<L;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
	for(int i=l;i<=mid;++i) a[i-l]=dp[i];
	for(int i=mid+1;i<L+l;++i)  a[i-l]=0;
	for(int i=0;i<L;++i) b[i]=fac[i];
	NTT(a,H);  NTT(b,H);
	for(int i=0;i<L;++i) a[i]=mul(a[i],b[i]);
	NTT(a,H,1);
	for(int i=mid+1;i<=r;++i) (dp[i]+=(mod-a[i-l]))%=mod; 
	CDQ(mid+1,r);
}
void pre(){
	w[0]=1;
	wn[0][18]=fast(G,(mod-1)/(1<<18));
	wn[1][18]=fast(wn[0][18],mod-2);
	for(int i=17;i>=0;--i)for(int k=0;k<=1;++k)wn[k][i]=mul(wn[k][i+1],wn[k][i+1]);
}
void Init(){
	fac[0]=1; for(int i=1;i<=100000;++i) fac[i]=mul(fac[i-1],i);
	dp[0]=0;  for(int i=1;i<=100000;++i) dp[i]=fac[i];
    pre();  CDQ(1,100000);
}
int solve(){
	n=read(); m=read(); int ans=1,flag=0;
	for(int i=1;i<=m;++i){
		int cnt=read(),maxx=0,minn=INF;
		for(int j=1;j<=cnt;++j){
			int x=read();
			cmax(maxx,x);
			cmin(minn,x);
		}
		if(maxx-minn+1!=cnt)flag=1;
		ans=mul(ans,dp[cnt]);
	}
	if(flag) return 0;
	else return ans;
}
int main(){
	//freopen(FILE".in","r",stdin);
	//freopen(FILE".out","w",stdout);
	Init(); for(int T=read();T;--T)printf("%d\n",solve());
	return 0;
}

例二

【HDU5322】Hope

题目大意：给出一个排列，每个点向它右面离它最近的且比它大的点连无向边，每个连通块的价值为这个连通块大小的平方，每个排列的价值为所有连通块之积。求n个数的所有排列的价值之和。

解析：考虑到放完前i个数再放i+1，不管i+1放在哪里，i+1之前的数形成了一个连通块，且对之后的没有影响。

于是可以列出方程

Dp[i]=sigmak(dp[i-k]k^2C(i-1,k-1)*(k-1)!)

Dp[i]/(i-1)!=sigmak(dp[i-k]/(i-k)!*k^2)

然后就可以CDQ分治做了

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 300100
using namespace std;
const int mod=998244353,G=10;
int n,w[MAXN],f[MAXN],fac[MAXN],A[MAXN],B[MAXN],R[MAXN],inv[MAXN],wn[2][20];
inline int read(){
	int x=0,f=1;  char ch=getchar();
	while(!isdigit(ch))  {if(ch=='-')  f=-1; ch=getchar();}
	while(isdigit(ch))  {x=x*10+ch-'0';  ch=getchar();}
	return x*f;
}
inline void swap(int &a,int &b){a=a^b,b=a^b,a=a^b;}
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline int sub(int a,int b){return (a-=b)<0?a+mod:a;}
inline int mul(int a,int b){return 1LL*a*b%mod;}
inline int fast(int a,int b){int ret(1);for(;b;b>>=1,a=mul(a,a))if(b&1)ret=mul(ret,a);return ret;}
void NTT(int *a,int H,int f=0){
	int L=(1<<H);
	for(int i=0;i<L;++i) if(i<R[i]) swap(a[i],a[R[i]]);
	for(int e=1;e<=H;++e){
		int len=(1<<e),l=(len>>1);
		for(int i=1;i<l;++i) w[i]=mul(w[i-1],wn[f][e]);
		for(int st=0;st<L;st+=len)for(int k=0;k<l;++k){
			int x=a[st+k],y=mul(w[k],a[st+k+l]);
			a[st+k]=add(x,y); a[st+k+l]=sub(x,y);
		}
	}
	if(f) for(int i=0;i<L;++i) a[i]=mul(a[i],fast(L,mod-2));
}
void CDQ(int l,int r){
	if(l==r) return;
	int mid=(l+r)>>1;  
	CDQ(l,mid);	
	int H=0,L=1; while(L<=r-l+1) L<<=1,++H;
	for(int i=0;i<L;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
	for(int i=l;i<=mid;++i) A[i-l]=mul(f[i],inv[i]);
	for(int i=mid+1;i<L+l;++i) A[i-l]=0;
	for(int i=0;i<L;++i) B[i]=mul(i,i);
	NTT(A,H);  NTT(B,H);
	for(int i=0;i<L;++i) A[i]=mul(A[i],B[i]);
	NTT(A,H,1);
	for(int i=mid+1;i<=r;++i) f[i]=add(f[i],mul(A[i-l],fac[i-1]));
	CDQ(mid+1,r);
}
void pre(){
	int N=100000;  w[0]=1;
	fac[0]=1; for(int i=1;i<=N;++i) fac[i]=mul(fac[i-1],i);
	inv[N]=fast(fac[N],mod-2); 
	for(int i=N-1;i>=0;--i) inv[i]=mul(inv[i+1],i+1);
	wn[0][18]=fast(G,(mod-1)/(1<<18));
	wn[1][18]=fast(wn[0][18],mod-2);
	for(int i=17;i>=0;--i)for(int k=0;k<=1;++k)wn[k][i]=mul(wn[k][i+1],wn[k][i+1]);
	f[0]=1;  CDQ(0,N);
}
int main(){
	//freopen(FILE".in","r",stdin);
	//freopen(FILE".out","w",stdout);
	pre(); while(~scanf("%d",&n)) printf("%d\n",f[n]); 
	return 0;
}

和数据结构套在一起

例一（分块FFT）

【bzoj3509】COUNTARI

考虑暴力做法，枚举中间元素，两边暴力FFT，时间复杂度为$O(n^2logn)$。

然而我们可以优化一些东西。

我们把序列分块，把有至少两个数在同一个块内的(i,j,k)暴力求出来,设块大小是B，n*B的，B大概要设到2000左右。然后对于三个点都在不同块的再暴力FFT就好了。

#include<bits/stdc++.h>
#define FILE "read"
#define MAXN 300010
#define N 30000
typedef long long ll;
const double pi=acos(-1);
ll n,block,a[MAXN],R[MAXN],pre[MAXN],suf[MAXN];
inline int read(){
    int x=0,f=1;  char ch=getchar();
    while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getchar();}
    while(isdigit(ch))  {x=x*10+ch-'0';  ch=getchar();}
    return x*f;
}
struct complex{
    double r,v;
    complex(double a=0,double b=0):r(a),v(b){}
    inline complex operator+(const complex &b){return complex(r+b.r,v+b.v);}
    inline complex operator-(const complex &b){return complex(r-b.r,v-b.v);}
    inline complex operator*(const complex &b){return complex(r*b.r-v*b.v,r*b.v+v*b.r);}
}A[MAXN],B[MAXN],C[MAXN],w[MAXN];
inline void swap(complex &a,complex &b){complex t(a);a=b;b=t;}
inline ll min(ll a,ll b){return a<b?a:b;}
inline ll max(ll a,ll b){return a<b?b:a;}
void FFT(complex *a,int L,int f){
    for(int i=0;i<L;++i) if(i<R[i]) swap(a[i],a[R[i]]);
    for(int len=2;len<=L;len<<=1){
        int l=(len>>1);  //complex wn(cos(pi/l),f*sin(pi/l));
        for(int i=1;i<l;++i)  w[i]=w[i-1]*wn;
        for(int st=0;st<L;st+=len)for(int k=0;k<l;++k){
            complex x=a[st+k],y=w[k]*a[st+k+l];
            a[st+k]=x+y;  a[st+k+l]=x-y;
        }
    }
    if(f==-1) for(int i=0;i<L;++i) a[i].r/=L;
}
int main(){
    //freopen(FILE".in","r",stdin);
    //freopen(FILE".out","w",stdout);
    n=read();  //block=sqrt/(n);  
    block=1448;
    ll ans=0;
    for(int i=1;i<=n;++i)  a[i]=read();
    for(int i=1;i<=n;++i)  suf[a[i]]++;
    for(int i=1;i<=n;i+=block){
        int l=i,r=min(i+block-1,n);
        for(int j=l;j<=r;++j)  suf[a[j]]--;
        for(int j=l;j<=r;++j){
            for(int k=j+1;k<=r;++k){
                if(a[j]*2-a[k]>=0) ans+=pre[a[j]*2-a[k]];
                if(a[k]*2-a[j]>=0) ans+=suf[a[k]*2-a[j]];
            }
            pre[a[j]]++;
        }
    }                
    int L=1,H=0;  w[0]=1;
    while(L<=N+N)L<<=1,++H;
    for(int i=0;i<L;++i) R[i]=(R[i>>1]>>1)|((i&1)<<(H-1));
    for(int i=1;i<=n;i+=block){
        int l=i,r=min(i+block-1,n);
        for(int j=0;j<=N;++j)    pre[j]=suf[j]=0;
        for(int j=1;j<l;++j)     pre[a[j]]++;
        for(int j=r+1;j<=n;++j)  suf[a[j]]++;
        for(int j=0;j<=N;++j)    A[j]=complex(pre[j],0);
        for(int j=0;j<=N;++j)    B[j]=complex(suf[j],0);
        for(int j=N+1;j<=L*2;++j)  A[j]=complex(0,0);
        for(int j=N+1;j<=L*2;++j)  B[j]=complex(0,0);
        FFT(A,L,1); FFT(B,L,1);
        for(int j=0;j<L;++j)  C[j]=A[j]*B[j];
        FFT(C,L,-1);
        for(int j=l;j<=r;++j)  ans+=(ll)(C[a[j]*2].r+0.5);
    }
    printf("%lld\n",ans);
    return 0;
}

本文标题:【算法笔记】多项式快速算法——FFT/NTT的应用

文章作者:chty

发布时间:2017年05月05日 - 20时34分

最后更新:2018年07月09日 - 20时19分

许可协议:本文为博主原创，未经许可不得转载。