既然选择了远方,便只顾风雨兼行|

H_W_Y

园龄:1年10个月粉丝:28关注:15

Poly 杂记(一)

比较杂,没有什么顺序。


板子

二话不说,先来一个板子。

Code
#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
using namespace __gnu_pbds;
using namespace std;
#define ll long long
#define ull unsigned long long

const int N=2.1e6+5,H=998244353;
int n,m,U[N],inv[N];
ull F[N],G[N],w[N];

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int dec(int a,int b){return a<b?a-b+H:a-b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

const int g=3,ig=qpow(g);

void init(){
  inv[1]=1;
  for(int i=2;i<N;i++) inv[i]=mul(inv[H%i],dec(0,H/i));
}

void Init(int n){//预处理蝴蝶变换
  for(int i=0;i<n;i++) U[i]=(U[i>>1]>>1)|((i&1)?n>>1:0);
}

ull* Cro(ull *F,ull *G,int n){//点乘
  for(int i=0;i<n;i++) F[i]=mul(F[i],G[i]);
  return F;
}

void Clr(ull *F,int l,int r){//[l,r) 清空
  for(int i=l;i<r;i++) F[i]=0;
}

ull* Times(ull *F,int n,int x){//多项式数乘
  for(int i=0;i<n;i++) F[i]=mul(F[i],x);
  return F;
}

ull* Dlt(ull *F,int n){//求导
  for(int i=1;i<n;i++) F[i-1]=mul(F[i],i);
  F[n-1]=0;return F;
}

ull* Sum(ull *F,int n){//导函数 -> 原函数
  for(int i=n;i>0;i--) F[i]=mul(F[i-1],inv[i]);
  F[0]=0;return F;
}

void Rev(ull *F,int l,int r){//翻转
  reverse(F+l,F+r);
}

void Cpy(ull *A,ull *B,int n){//赋值
  memcpy(A,B,n<<3);
}

void NTT(ull *p,int len,int op){//DFT/IDFT
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  ull u,tmp;w[0]=1;
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=qpow(op==1?g:ig,(H-1)/k);
	for(int i=1;i<=lst;i++) w[i]=mul(w[i-1],u);
    for(int i=0;i<len;i+=k)
	  for(int l=i;l<i+lst;l++)
	    tmp=w[l-i]*p[l+lst]%H,p[l+lst]=p[l]-tmp+H,p[l]+=tmp;
    if(k==(1<<16)) for(int i=0;i<len;i++) p[i]%=H;
  }
  for(int i=0;i<len;i++) p[i]%=H;
}

ull* Mul(ull *F,ull *G,int n,int m){//多项式乘法
  int P=1,I;
  for(P=1;P<n+m;P<<=1);
  Clr(F,n,P),Clr(G,m,P);
  I=qpow(P),Init(P);
  NTT(F,P,1),NTT(G,P,1),Cro(F,G,P),NTT(F,P,-1);
  Times(F,P,I);
  return F;
}

ull* Inv(ull *F,int n){//多项式乘法逆
  int P=1,I;
  for(P=1;P<n;P<<=1);
  ull A[N],B[N],C[N];
  Clr(A,0,P<<2),Clr(B,0,P<<2),Clr(C,0,P<<2);

  A[0]=qpow(F[0]);
  for(int k=2;k<=P;k<<=1){
	Cpy(B,F,k),Cpy(C,A,k>>1);
	Clr(A,k>>1,k),Init(k),I=qpow(k);

	NTT(A,k,1),NTT(B,k,1),Cro(B,A,k),NTT(B,k,-1);
	Times(B,k,I);

	Clr(B,0,k>>1),NTT(B,k,1),Cro(B,A,k),NTT(B,k,-1);
	Times(B,k,I);

	for(int i=0;i<(k>>1);i++) A[i]=C[i];
	for(int i=(k>>1);i<k;i++) A[i]=dec(0,B[i]);
  }
  for(int i=0;i<n;i++) F[i]=A[i];
  return F;
}

ull* Ln(ull *F,int n){//多项式 ln(F[0]=1)
  ull A[N];
  Clr(A,0,n<<1),Cpy(A,F,n);
  Inv(A,n),Dlt(F,n),Mul(F,A,n,n);
  Sum(F,n-1),Clr(F,n,n<<1);
  return F;
}

ull* Exp(ull *F,int n){//多项式 exp(F[0]=0)
  int P=1;
  for(P=1;P<n;P<<=1);
  ull A[N],B[N];
  A[0]=1,Clr(A,1,P<<1),Clr(B,0,P<<1);
  for(int k=2;k<=P;k<<=1){
	Cpy(B,A,k>>1),Clr(B,k>>1,k<<1),Ln(B,k);
	for(int i=0;i<k;i++) B[i]=(F[i]-B[i]+H)%H;
	B[0]=(B[0]+1)%H,Mul(A,B,k,k),Clr(A,k,k<<1);
  }
  for(int i=0;i<n;i++) F[i]=A[i];
  return F;
}

ull* Qpow(ull *F,int n,int k){//多项式快速幂(F[0]=1)
  Ln(F,n),Times(F,n,k),Exp(F,n);
  return F;
}

ull* Qpow_plus(ull *F,int n,int k1,int k2){//多项式快速幂加强版(F[0] 无限制,在外面要把全是 0 的情况判掉)
  ull A[N],id=0;Clr(A,0,n<<1);
  for(int i=0;i<n;i++) if(F[i]!=0){
    int I=qpow(F[i]);
	for(int j=i;j<n;j++)
	  A[j-i]=mul(F[j],I);
	id=i;
	break;
  }
  Qpow(A,n-id,k1);
  int nw=qpow(F[id],k2);
  Clr(F,0,n),id*=(ll)k1;
  for(ull i=id;i<(ull)n;i++) F[i]=mul(A[i-id],nw);
  return F;
}

ull* Sqrt_1(ull *F,int n,int k){//多项式 k 次方根(F[0]=1),Ln + Exp 实现,常数较大
  Ln(F,n),Times(F,n,inv[k]),Exp(F,n);
  return F;
}

int Sqrt_BSGS(int a,int b){
  int A=1,B=sqrt(H)+1;a%=H,b%=H;
  if(H==1||b==1) return 0;
  if(!a) return !b?1:-1;
  gp_hash_table<int,int> mp;
  for(int i=1;i<=B;i++) mp[mul(A=mul(A,a),b)]=i;
  for(int i=1,AA=A;i<=B;i++,AA=mul(AA,A))
    if(mp.find(AA)!=mp.end()) return i*B-mp[AA];
  return -1;
}

ull Sqrt_H(int x,int k=2){
  int res=qpow(g,Sqrt_BSGS(g,x)/k);
  return min(res,H-res);
}

ull* Sqrt(ull *F,int n,int k=2){//多项式 k 次方根
  int P=1,I=qpow(2);
  for(P=1;P<n;P<<=1);
  ull A[N],B[N],C[N];
  Clr(A,0,P<<1),Clr(B,0,P<<1),Clr(C,0,P<<1);
  A[0]=Sqrt_H(F[0],k);

  for(int k=2;k<=P;k<<=1){
    Cpy(B,A,k>>1),Clr(B,k>>1,k),Inv(B,k),Clr(B,k,k<<1);
	Cpy(C,F,k),Mul(B,C,k,k);
	for(int i=0;i<(k<<1);i++) B[i]=(B[i]+A[i])%H;
	for(int i=0;i<k;i++) A[i]=mul(B[i],I);
  }
  for(int i=0;i<n;i++) F[i]=A[i];
  return F;
}

ull* Mof(ull *F,ull *G,int n,int m){//多项式带余除法
  ull A[N],B[N];
  Clr(A,0,n+m),Clr(B,0,n+m);
  int P=n-m+1;
  for(int i=m;i<n;i++) G[i]=0;
  Rev(F,0,n),Cpy(A,F,P),Rev(F,0,n);
  Rev(G,0,m),Cpy(B,G,P),Rev(G,0,m);
  Inv(B,P),Mul(B,A,P,P),Rev(B,0,P),Clr(B,P,P<<1),Cpy(A,B,P);
  Mul(B,G,n,n),Clr(B,n,n<<1);
  for(int i=0;i<m-1;i++) G[i]=(F[i]-B[i]+H)%H;
  for(int i=0;i<P;i++) F[i]=A[i];
  G[m-1]=0,Clr(F,P,P<<1);
  return F;
}

ull* Sin(ull *F,int n){//多项式 Sin 函数
  ull I=Sqrt_H(H-1),A[N];Clr(A,0,n<<1);
  Times(F,n,I),Exp(F,n);
  Cpy(A,F,n),Inv(A,n),I=qpow(I<<1);
  for(int i=0;i<n;i++) F[i]=I*(F[i]-A[i]+H)%H;
  return F;
}

ull* Cos(ull *F,int n){//多项式 Cos 函数
  ull I=Sqrt_H(H-1),A[N];Clr(A,0,n<<1);
  Times(F,n,I),Exp(F,n);
  Cpy(A,F,n),Inv(A,n),I=qpow(2);
  for(int i=0;i<n;i++) F[i]=I*(F[i]+A[i])%H;
  return F;
}

#define mid ((l+r)>>1)
#define lc p<<1
#define rc p<<1|1
#define lson l,mid,lc
#define rson mid+1,r,rc

namespace Evaluation{//多项式多点求值
  ull *Q[N],tr[N<<5],*np(tr),A[N],B[N];
  int len[N];

  void build(int l,int r,int p,ull *a){
    Q[p]=np,np+=r-l+2;len[p]=r-l+2;
	if(l==r) return Q[p][1]=1,Q[p][0]=(H-a[l])%H,void();
	build(lson,a),build(rson,a);

	Cpy(A,Q[lc],len[lc]),Cpy(B,Q[rc],len[rc]);
	Mul(A,B,len[lc],len[rc]);
    for(int i=0;i<len[p];i++) Q[p][i]=A[i]; 
  }

  void solve(int l,int r,int p,ull *G){
    if(l==r) return G[l]=Q[p][0],void();
	Cpy(A,Q[p],(len[p]-1)),Cpy(B,Q[lc],len[lc]);
	Mof(A,B,len[p]-1,len[lc]);
	Cpy(Q[lc],B,len[lc]);

	Cpy(A,Q[p],len[p]),Cpy(B,Q[rc],len[rc]);
	Mof(A,B,len[p]-1,len[rc]);
	Cpy(Q[rc],B,len[rc]);
	solve(lson,G),solve(rson,G);
  }

  ull* Multipoint(ull *F,ull *a,int n,int m,ull *G){
    build(0,m-1,1,a);
	if(n>=len[1]){
	  Cpy(A,Q[1],len[1]);
	  Mof(F,A,n,len[1]);
	  Cpy(Q[1],A,len[1]);
	}else Cpy(Q[1],F,len[1]);
    solve(0,m-1,1,G);
	return G;
  }
}

namespace Interpolation{//多项式快速插值
  ull *Q[N],*P[N],tr[N<<6],*np(tr),A[N],B[N],R[N],C[N],D[N];
  int len[N];

  void build(int l,int r,int p,ull *a){
    Q[p]=np,np+=r-l+2,P[p]=np,np+=r-l+2;len[p]=r-l+2;
	if(l==r) return Q[p][1]=1,Q[p][0]=(H-a[l])%H,void();
	build(lson,a),build(rson,a);

	Cpy(A,Q[lc],len[lc]),Cpy(B,Q[rc],len[rc]);
	Mul(A,B,len[lc],len[rc]);
    for(int i=0;i<len[p];i++) Q[p][i]=A[i]; 
  }

  void solve(int l,int r,int p){
    if(l==r) return P[p][0]=R[l],P[p][1]=0,void();
	solve(lson),solve(rson);
	Cpy(A,Q[lc],len[lc]),Cpy(B,Q[rc],len[rc]);
	Cpy(C,P[lc],len[lc]),Cpy(D,P[rc],len[rc]);
    Mul(A,D,len[lc],len[rc]-1),Mul(B,C,len[rc],len[lc]-1);
	for(int i=0;i<len[p];i++) P[p][i]=(A[i]+B[i])%H;
  }

  void Lagrange(ull *a,ull *b,int n,ull *G){
	build(0,n-1,1,a);
	Dlt(Q[1],len[1]);
    Evaluation::Multipoint(Q[1],a,n,n,R);
    for(int i=0;i<n;i++) R[i]=mul(qpow(R[i]),b[i]);
    solve(0,n-1,1);
	Cpy(G,P[1],n);
  }
}

namespace Partition_FFT{//分治 FFT
  ull A[N],B[N];
  
  void solve(int l,int r,ull *F,ull *G){
	if(l==r) return;
	solve(l,mid,F,G);
	Cpy(A,G+l,mid-l+1),Cpy(B,F,r-l+1);
	Mul(A,B,mid-l+1,r-l+1);
	for(int i=mid+1;i<=r;i++) G[i]=(G[i]+A[i-l])%H;
	solve(mid+1,r,F,G);
  }

  void Partition(ull *F,int n,ull *G){
    Clr(G,0,n),G[0]=1;
	solve(0,n-1,F,G);
  }
}

namespace FWT{
  ull A[N],B[N];
  
  void Or(ull *p,int len,int op){
	for(int k=2,lst=1;k<=len;lst=k,k<<=1)
	  for(int i=0;i<len;i+=k)
	    for(int l=i;l<i+lst;l++)
		  (p[l+lst]+=op*p[l]+H)%=H;
  }

  void And(ull *p,int len,int op){
	for(int k=2,lst=1;k<=len;lst=k,k<<=1)
	  for(int i=0;i<len;i+=k)
	    for(int l=i;l<i+lst;l++)
		  (p[l]+=op*p[l+lst]+H)%=H;
  }

  void Xor(ull *p,int len,int op){
	for(int k=2,lst=1;k<=len;lst=k,k<<=1)
	  for(int i=0;i<len;i+=k)
	    for(int l=i;l<i+lst;l++){
		  (p[l]+=p[l+lst])%=H,p[l+lst]=(p[l]-2ll*p[l+lst]%H+H)%H;
		  p[l]=1ll*p[l]*op%H,p[l+lst]=1ll*p[l+lst]*op%H;
		}
  }

  void Mul_or(ull *F,ull *G,int n){
	Cpy(A,F,n),Cpy(B,G,n);
	Or(A,n,1),Or(B,n,1);
	Cro(A,B,n);
	Or(A,n,-1);
	for(int i=0;i<n;i++) cout<<A[i]<<' ';
	cout<<'\n';
  }

  void Mul_and(ull *F,ull *G,int n){
	Cpy(A,F,n),Cpy(B,G,n);
	And(A,n,1),And(B,n,1);
	Cro(A,B,n);
	And(A,n,-1);
	for(int i=0;i<n;i++) cout<<A[i]<<' ';
	cout<<'\n';
  }

  void Mul_xor(ull *F,ull *G,int n){
	Cpy(A,F,n),Cpy(B,G,n);
	Xor(A,n,1),Xor(B,n,1);
	Cro(A,B,n);
	Xor(A,n,qpow(2));
	for(int i=0;i<n;i++) cout<<A[i]<<' ';
	cout<<'\n';
  }

}

int main(){
  init();
  //...
  return 0;
}

一下子大小变成了 9k。

以下的例题有些是没有放板子的,注意一下 清空 就可以了。但是每一遍都是重新写的。


有用的式子

ln(1x)=i=1xii

证明两边同时求导即可。


下面例题中有些使用 vector 实现的多项式操作,具体存在于 CF755G 和 UOJ182。


*MTT

感觉还是需要单独提一嘴,后面单位根反演部分几乎是道道题要 MTT,真是令人难绷。/qd


其核心是把多项式系数 Fi 拆成 Aik+Bi 的形式,从而解决 值域过大精度过高 的问题。

例如在 P4245 【模板】任意模数多项式乘法 中,乘出来的值高达 109×109×105=1023,我们尝试取 k=215,那么就可以把 A,B 的最终值域优化到 1014,非常优秀。

那么这道题就变成了求

(A0k+A1)(B0k+B1)=A0B0k2+(A0B1+A1B0)k+A1B1

直接做可以做到用 7 次 FFT,即四次 DFT,三次 IDFT,可以跑出不错的成绩。Record

Code
#include <bits/stdc++.h>
using namespace std;
#define ld long double
#define ll long long
#define cp complex<ld>

const int N=2.1e6+3,Bl=(1<<15);
const ld PI=acos(-1);
int n,m,P,U[N];
cp A0[N],A1[N],B0[N],B1[N],A[N],B[N],C[N];
ll ans=0,H;

void FFT(cp *p,int len,int op){
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  cp nw=cp(1,0),tmp,u;
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
    u=cp(cos(2*PI/k),op*sin(2*PI/k));
	for(int i=0;i<len;i+=k){
	  nw=cp(1,0);
	  for(int l=i;l<i+lst;l++,nw=nw*u)
	    tmp=nw*p[l+lst],p[l+lst]=p[l]-tmp,p[l]=p[l]+tmp;
	}
  }
}

int main(){
  /*2024.4.29 H_W_Y P4245 【模板】任意模数多项式乘法 MTT*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>m>>H;

  for(P=1;P<=n+m;P<<=1);
  for(int i=0;i<P;i++) U[i]=(U[i>>1]>>1)|((i&1)?P>>1:0);

  for(int i=0,x;i<=n;i++) cin>>x,A0[i]=x/Bl,A1[i]=x%Bl;
  for(int i=0,x;i<=m;i++) cin>>x,B0[i]=x/Bl,B1[i]=x%Bl;
  FFT(A0,P,1),FFT(A1,P,1),FFT(B0,P,1),FFT(B1,P,1);

  for(int i=0;i<P;i++)
    A[i]=A0[i]*B0[i],B[i]=A0[i]*B1[i]+A1[i]*B0[i],C[i]=A1[i]*B1[i];
  FFT(A,P,-1),FFT(B,P,-1),FFT(C,P,-1);

  for(int i=0;i<=n+m;i++){
    ll a,b,c,ans;
	a=(ll)(A[i].real()/P+0.5),b=(ll)(B[i].real()/P+0.5),c=(ll)(C[i].real()/P+0.5);
    ans=(1ll*a%H*Bl%H*Bl%H+1ll*b%H*Bl%H+c%H)%H;
	cout<<ans<<' ';
  }
  return 0;
}

对这个东西我们还可以进行进一步的优化,也就是类似于 FFT 中的三步变两步,将其转化成复数多项式。

也就是我们可以设

P=A0+A1iQ=B0+B1iP=A0A1i

那么就可以得到

PQ=A0B0A1B1+(A0B1+A1B0)iPQ=A0B0+A1B1+(A0B1A1B0)i

于是我们需要的东西都可以从上面两个东西做加减出来了,容易发现这样只需要 5 次 FFT。

果不其然会比上面快不少。Record

Code
#include <bits/stdc++.h>
using namespace std;
#define ld long double
#define cp complex<ld>
#define ll long long

const int N=2.1e6+3,Bl=(1<<15);
const ld PI=acos(-1);
int n,m,P,U[N];
ll H;
cp F[N],G[N],A[N],B[N],C[N];

void FFT(cp *p,int len,int op){
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  cp nw,tmp,u;
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=cp(cos(2*PI/k),op*sin(2*PI/k));
	for(int i=0;i<len;i+=k){
	  nw=cp(1,0);
	  for(int l=i;l<i+lst;l++,nw=nw*u)
	    tmp=nw*p[l+lst],p[l+lst]=p[l]-tmp,p[l]=p[l]+tmp;
	}
  } 
}

int main(){
  /*2024.4.29 H_W_Y P4245 【模板】任意模数多项式乘法 MTT-5*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>m>>H;

  for(P=1;P<=n+m;P<<=1);
  for(int i=0;i<P;i++) U[i]=(U[i>>1]>>1)|((i&1)?P>>1:0);

  for(int i=0,x;i<=n;i++) cin>>x,A[i]+=cp(x/Bl,x%Bl),C[i]+=cp(x/Bl,-x%Bl);
  for(int i=0,x;i<=m;i++) cin>>x,B[i]+=cp(x/Bl,x%Bl);

  FFT(A,P,1),FFT(B,P,1),FFT(C,P,1);
  for(int i=0;i<P;i++) F[i]=A[i]*B[i],G[i]=C[i]*B[i];
  FFT(F,P,-1),FFT(G,P,-1);

  for(int i=0;i<=n+m;i++){
	ll a,b,c,ans;
	a=(ll)((F[i].real()+G[i].real())/2/P+0.5),b=(ll)(F[i].imag()/P+0.5),c=(ll)((G[i].real()-F[i].real())/2/P+0.5);
	ans=(1ll*a%H*Bl%H*Bl%H+1ll*b%H*Bl%H+c%H)%H;
    cout<<ans<<' ';
  }
  return 0;
}

P3338 [ZJOI2014] 力 - FFT

给出 n 个数 q1,q2,qn,定义

Fj = i=1j1qi×qj(ij)2  i=j+1nqi×qj(ij)2,Ei = Fiqi

1in,求 Ei 的值。

1n1050<qi<109

典题先来一个。


首先推一下式子

Ej=i<jqi(ij)2i>jqi(ij)2

我们设

dk={1k2k>00k=01k2k<0

于是我们有

Ej=i=0nqidji

但是这不是标准卷积形式,所以我们考虑 对称(镜像)。

Ej+n=Ej,Dk+n=Dk,于是就有了以下式子

Ej+n=i+k=j+nqiDk

这就是一个卷积形式了,直接 FFT 就可以了。

注意精度问题,我们需要给 D 乘上 1011代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ld long double
#define cp complex<ld>

const int N=2.1e6+5;
const ld PI=acos(-1),B=1e11;
int n,m,U[N];
cp F[N];
ld a[N],b[N];

void FFT(cp *p,int len,int op){
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  cp u,tmp,nw=cp(1,0);
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=cp(cos(2*PI/k),op*sin(2*PI/k));
	for(int i=0;i<len;i+=k){
	  nw=cp(1,0);
	  for(int l=i;l<i+lst;l++,nw=nw*u)
	    tmp=nw*p[l+lst],p[l+lst]=p[l]-tmp,p[l]+=tmp;
	}
  }
}

void Mul(ld *a,ld *b,int n,int m){
  int P=1;
  for(P=1;P<n+m;P<<=1);
  for(int i=0;i<P;i++) U[i]=(U[i>>1]>>1)|((i&1)?P>>1:0);

  for(int i=0;i<P;i++) F[i]=cp(a[i],b[i]);
  FFT(F,P,1);
  for(int i=0;i<P;i++) F[i]=F[i]*F[i];
  FFT(F,P,-1);
  for(int i=0;i<P;i++) a[i]=F[i].imag()/P/2.0;
}

int main(){
  /*2024.7.10 H_W_Y P3338 [ZJOI2014] 力 FFT*/
  scanf("%d",&n);
  for(int i=0;i<n;i++) scanf("%Lf",&a[i]);
  for(int i=0;i<=2*n;i++){
    if(i-n>0) b[i]=(ld)1.0/(i-n)/(i-n);
	else if(i-n<0) b[i]=(ld)-1.0/(i-n)/(i-n);
	b[i]*=B;
  }
  Mul(a,b,n,2*n);
  for(int i=n;i<(n<<1);i++) printf("%.3Lf\n",a[i]/B);
  return 0;
}

UVA12633 Super Rooks on Chessboard - 容斥

R×C 的棋盘,有一些”超级车“,问:有多少个格子未被攻击。

一辆车可以攻击其所在行,列和对角线。

R,C,N105

小容斥一手。


考虑把问题转化成求有多少个格子被至少一辆车攻击。

于是我们设

{R:C:D:线

那么其实答案是

|RCD|=|R|+|C|+|D||RC||RD||CD|+|RCD|

前面的 |R|,|C|,|D|,|RC|,|RD|,|CD| 都是可以在 O(1)/O(n) 的时间复杂度内求出来的,而且是简单的。

于是其实只用讨论如何求 |RCD| 了。


不妨设

{r1,r2,,rnc1,c2,,cmd1,d2,,dl线

合理地对对角线编号,问题转化成求有多少 (i,j,k) 使得 cjri=dk,移一下项就可以得到

ri+dk=cj

这种和式不就是 多项式卷积 形式?!

于是我们令

R(x)=xri,D(x)=xdk

那么答案就是

|RCD|=j[xcj]R(x)D(x)

这是可以直接 FFT 之后线性求出来的,结束!

时间复杂度 O(nlogn)代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ld double
#define ll long long
#define cp complex<ld>

const int N=1e6+5;
const ld PI=acos(-1);
int R,C,num,n,m,l,r[N],c[N],d[N],U[N];
ll ans=0,a[N],b[N];
bool vr[N],vc[N],vd[N];
cp F[N];

void FFT(cp *p,int len,int op){}
void Mul(ll *a,ll *b,int n,int m){}

void SOLVE(int id){
  cin>>R>>C>>num;
  memset(vr,0,(R+1)<<1),memset(vc,0,(C+1)<<1),memset(vd,0,(R+C+1)<<1);
  memset(a,0,(R+C+1)<<3),memset(b,0,(R+C+1)<<3);
  for(int i=0,x,y;i<num;i++) cin>>x>>y,vr[x]=1,vc[y]=1,vd[y-x+R]=1;

  n=m=l=0;
  for(int i=1;i<=R;i++) if(vr[i]) r[++n]=i;
  for(int i=1;i<=C;i++) if(vc[i]) c[++m]=i;
  for(int i=1;i<=R+C;i++) if(vd[i]) d[++l]=i;  
  
  ans=1ll*n*C+1ll*m*R-1ll*n*m;
  for(int i=1;i<=l;i++){
	if(d[i]-R>=0) ans+=min(R,C-d[i]+R);
	else ans+=min(C,d[i]);
  }

  for(int i=1,it1=l+1,it2=l;i<=n;i++){
	while(it1>1&&d[it1-1]>=1-r[i]+R) --it1;
    while(it2>=1&&d[it2]>C-r[i]+R) --it2;
    ans-=(it2-it1+1);
  }

  for(int i=1,it1=1,it2=0;i<=m;i++){
	while(it1<=l&&d[it1]<c[i]) ++it1;
	while(it2<l&&d[it2+1]<=c[i]-1+R) ++it2;
	ans-=(it2-it1+1);
  }

  for(int i=1;i<=n;i++) a[r[i]]=1;
  for(int i=1;i<=l;i++) b[d[i]]=1;

  Mul(a,b,R+1,R+C);

  for(int i=1;i<=m;i++) ans+=a[c[i]+R];
  cout<<"Case "<<id<<": "<<(1ll*R*C-ans)<<'\n';
}

int main(){
  /*2024.7.10 H_W_Y UVA12633 Super Rooks on Chessboard FFT*/
  ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
  int _;cin>>_;
  for(int i=1;i<=_;i++) SOLVE(i);
  return 0;
}

SP8372 TSUM - Triple Sums - 小容斥

N个整数A1,A2,,AN,对所有的S,求满足$A_i+A_j+A_k = S i < j < k(i,j,k)$对数.

N40000,Ai20000


首先来考虑如果没有 i,j,k 的限制怎么办?

容易发现其实就是

[xS](i=1NxAi)3

但是这里面有哪些重复的部分呢?

这个东西其实是

(i=1NxAi)3=i=1Nx3Ai+3ijx2Ai+Aj+6i<j<kxAi+Aj+Ak

最后的部分是我们想要求的。

而对于第一部分,我们可以 O(n) 求出;而对于第二项,我们就考虑再用一次生成函数,也就是用

(i=1Nx2Ai)(i=1NxAi)

再来分讨一下它是什么?

容易发现其实就是

(i=1Nx2Ai)(i=1NxAi)=i=1Nx3Ai+ijx2Ai+Aj

这样通过两次 FFT 我们就得到了想要得到的东西了。

于是就做完了,时间复杂度 O(VlogV)代码

Code
const int N=2e6+5,S=2e4+1;
const ld PI=acos(-1);
int n,m,A[N],mx,U[N];
cp F[N];
ll a[N],b[N],c[N];

void FFT(cp *p,int len,int op){}
void Mul(ll *a,ll *b,int n,int m){}

int main(){
  /*2024.7.11 H_W_Y SP8372 TSUM - Triple Sums FFT*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n;
  for(int i=0;i<n;i++) cin>>A[i],A[i]+=S;
  for(int i=0;i<n;i++) a[A[i]]++,b[2*A[i]]++,c[A[i]]++;

  Mul(b,a,S<<2,S<<1);
  Mul(c,a,S<<1,S<<1),Mul(c,a,S<<2,S<<1);

  for(int i=0;i<n;i++) a[A[i]]--,a[3*A[i]]++;

  for(int i=0;i<(S<<3);i++) a[i]=(c[i]-3ll*b[i]+2ll*a[i])/6ll;

  for(int i=0;i<(S<<3);i++) if(a[i])
	cout<<(i-3*S)<<" : "<<a[i]<<'\n';
  return 0;
}

 AT icpc2013spring F.Point Distance - 二维转一维

给定一张 N×N 的网格图。一些点会落在网格图上,记 Cxy 为坐标是 (x,y) 的点的数量 (0Cxy9)。现要求对于所有 成对的点之间的欧几里得距离,统计:

  • 全部距离依次平方后,排序并去重后的结果 d1,d2,,dm
  • 上述序列去重前 di 重复出现的次数 ci
  • 所有距离的平均距离 Dave

1N10240Cxy9

应该是紫吧?


先写一个式子出来

Dx,y=i=0N1j=0N1Ci,jCi+x,j+y

这是二维的式子,非常不好处理,所以我们考虑一种映射关系

ϕ(x,y)=x×2N+y

同时令 M=2N2,那么就可以得到

Dϕ(x,y)=i=0N1j=0N1Cϕ(i,j)Cϕ(i,j)+ϕ(x,y)

于是有

Dz=u=0MCuCu+z

这个东西很像卷积,但是又不是标准的卷积形式。

按照之前的套路,我们直接对它 翻转,也就是令

{Ck=CMkDk=DMk

于是就得到了

Dk=i=0kCiCki

这就可以用卷积直接做了,上 FFT 就完事,不卡精度不卡空间(表示好评)。

注意特殊处理一下 dis=0 的情况就可以了。

时间复杂度 O(MlogM)代码

Code
const int N=8.4e6+5;
const ld PI=acos(-1);
int n,m,U[N],num;
ll a[N],b[N],d[N],cnt;
cp F[N];
ld ans=0;

void FFT(cp *p,int len,int op){}
void Mul(ll *a,ll *b,int n,int m){}

int id(int x,int y){return x*2*n+y;}

int main(){
  /*2024.7.11 H_W_Y AT_icpc2013spring_f Point Distance FFT*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n;m=2*n*n;
  for(int i=0;i<n;i++)
    for(int j=0;j<n;j++)
	  cin>>a[id(i,j)],cnt+=a[id(i,j)];
  for(int i=0;i<=m;i++) b[m-i]=a[i];

  Mul(a,b,m+1,m+1);

  for(int i=0;i<=m;i++){
	int x=(m-i)/(2*n),y=(m-i)%(2*n);
	if(y>n) y=2*n-y,++x;
	// cerr<<x<<" "<<y<<endl;
	d[x*x+y*y]+=a[i];
  }
  d[0]-=cnt,d[0]/=2;
  cnt=0;
  for(int i=0;i<=m;i++) cnt+=d[i];
  for(int i=0;i<=m;i++) ans+=(ld)sqrt((ld)i)*d[i]/cnt;
  //ans/=(ld)cnt;

  cout<<fixed<<setprecision(10)<<ans<<'\n';
  num=0;
  for(int i=0;i<=m;i++){
	if(d[i]) cout<<i<<' '<<d[i]<<'\n',++num;
	if(num==10000) break;
  }
  return 0;
}

UVA12298 Super Poker II - FFT

有一副超级扑克,包含无数张牌。对于每个正合数 p,恰好有 4 张牌:黑桃 p,红桃 p,梅花 p 和方块 p(分别用 pSpHpCpD 表示)。没有其它类型的牌。

给定一个整数 n,从 4 种花色中各选一张牌,问有多少种组合可以使得点数之和 =n,例如,n=24 时,有一种组合方法是 4S+6H+4C+10D,如下图所示。

390770

不巧的是,有些牌已经丢失,并且为了让题目更毒瘤,我们还会提供两个正整数 ab,你的任务是按顺序输出 n=a,n=a+1,n=a+2,...,n=b 时的答案。

b50000,c10000

唐唐板子题。


直接把四个多项式乘起来就可以了。

注意每组数据之后要输出一个 \n,否则会喜提 presentation error。代码

Code
const int N=4.2e6+5;
const ld PI=acos(-1);
int n,m,U[N],l,r,cnt=0,p[N];
cp F[N];
ll a[N],b[N],c[N],d[N];
bool vis[N];
string s;

void FFT(cp *p,int len,int op){}
void Mul(ll *a,ll *b,int n,int m){}

void init(){
  for(int i=2;i<=50000;i++){
	if(!vis[i]) p[++cnt]=i;
	for(int j=1;j<=cnt&&i*p[j]<=50000;j++){
	  vis[p[j]*i]=true;
	  if(i%p[j]==0) break;
	}
  }
}

int gnum(string s){
  if((int)s.size()>6) return 0;
  int x=0;
  for(int i=0;i<(int)s.size()-1;i++) x=x*10+(s[i]-'0');
  return x;
}

void SOLVE(){
  cin>>l>>r>>n;m=r+1;
  if(!l&&!r&&!n) exit(0);

  memset(a,0,sizeof(a)),memset(b,0,sizeof(b)),memset(c,0,sizeof(c)),memset(d,0,sizeof(d));

  for(int i=2;i<=r;i++){
	if(!vis[i]) a[i]=b[i]=c[i]=d[i]=0;
    else a[i]=b[i]=c[i]=d[i]=1;
  }

  while(n--){
    cin>>s;char ch=s[(int)s.size()-1];
	if(ch=='S') a[gnum(s)]=0;
	else if(ch=='H') b[gnum(s)]=0;
    else if(ch=='C') c[gnum(s)]=0;
	else d[gnum(s)]=0;
  }

  Mul(a,b,m,m),Mul(c,d,m,m);
  Mul(a,c,m<<1,m<<1);

  for(int i=l;i<=r;i++) cout<<a[i]<<'\n';
  cout<<'\n';
}

int main(){
  /*2024.7.11 H_W_Y UVA12298 Super Poker II FFT*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  init();
  while(1) SOLVE();
  return 0;
}

CF755G PolandBall and Many Other Balls - 多项式优化 dp

有一排 n 个球,定义一个组可以只包含一个球或者包含两个相邻的球。现在一个球只能分到一个组中,求从这些球中取出 k 组的方案数。

n109k<215


首先,先写一个 dp,即设 fi,j 表示前 i 个数,分成 j 组方案数。

于是我们有

fi,j=fi1,j+fi1,j1+fi2,j1

由于需要 1k,再加上一些神秘原因,我们考虑找到 fi 的生成函数 Fi(x),那么就有

Fi(x)=(x+1)Fi1(x)+xFi2(x)

我的天,矩阵,快速,幂?

于是构造矩阵

A=[x+11x0]

那么答案其实就是 An[0][0]

做完了,直接暴力 NTT 即可,常数大,时间复杂度 O(klognlogk)

笔者本以为它过不了——但是它过了!需要开 C++20。代码

Code
void init(){
  for(P=1;P<(K<<1);P<<=1);
  I=qpow(P);
  for(int i=0;i<P;i++) U[i]=(U[i>>1]>>1)|((i&1)?P>>1:0);
}

poly operator *(poly a,poly b){
  poly A=a,B=b;
  A.resize(P),B.resize(P);

  NTT(A,P,1),NTT(B,P,1);
  for(int i=0;i<P;i++) A[i]=A[i]*B[i]%H;
  NTT(A,P,-1);

  A.resize(K);
  for(int i=0;i<K;i++) A[i]=1ll*A[i]*I%H;
  return A;
}

poly operator +(poly a,poly b){
  poly A(K,0);
  for(int i=0;i<K;i++) A[i]=(a[i]+b[i])%H;
  return A;
}

struct Mat{
  poly g[2][2];

  void init(){
	g[0][0]={1},g[0][0].resize(K);
	g[1][1]={1},g[1][1].resize(K);
	g[0][1].resize(K),g[1][0].resize(K);
  }
}a;

Mat operator *(Mat a,Mat b){
  Mat c;
  c.g[0][0]=(a.g[0][0]*b.g[0][0])+(a.g[0][1]*b.g[1][0]);
  c.g[0][1]=(a.g[0][0]*b.g[0][1])+(a.g[0][1]*b.g[1][1]);
  c.g[1][0]=(a.g[1][0]*b.g[0][0])+(a.g[1][1]*b.g[1][0]);
  c.g[1][1]=(a.g[1][0]*b.g[0][1])+(a.g[1][1]*b.g[1][1]);
  return c;
}

Mat Qpow(Mat a,int b){
  Mat c;c.init();
  while(b){if(b&1) c=c*a;a=a*a,b>>=1;}
  return c;
}

int main(){
  /*2024.7.11 H_W_Y CF755G PolandBall and Many Other Balls 不知道在干什么*/
  ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
  cin>>n>>K;++K;init();
  a.g[0][0]={1,1},a.g[0][0].resize(K);
  a.g[0][1]={1},a.g[0][1].resize(K);
  a.g[1][0]={0,1},a.g[1][0].resize(K);
  a.g[1][1]={0},a.g[1][1].resize(K);
  a=Qpow(a,n);
  for(int i=1;i<K;i++) cout<<a.g[0][0][i]<<' ';
  cout<<'\n';
}

注意到矩阵乘法部分是可以优化的,我们发现可以先把 8 次 DFT 做完之后,再点值相乘,最后做 4 次 IDFT。

这样可以把常数减小一半。代码

Code
Mat operator *(Mat a,Mat b){
  Mat c;
  for(int i=0;i<2;i++) for(int j=0;j<2;j++) c.g[i][j].resize(P),a.g[i][j].resize(P),b.g[i][j].resize(P),NTT(a.g[i][j],P,1),NTT(b.g[i][j],P,1);

  for(int i=0;i<P;i++){
    c.g[0][0][i]=(a.g[0][0][i]*b.g[0][0][i]+a.g[0][1][i]*b.g[1][0][i])%H;
    c.g[0][1][i]=(a.g[0][0][i]*b.g[0][1][i]+a.g[0][1][i]*b.g[1][1][i])%H;
    c.g[1][0][i]=(a.g[1][0][i]*b.g[0][0][i]+a.g[1][1][i]*b.g[1][0][i])%H;
    c.g[1][1][i]=(a.g[1][0][i]*b.g[0][1][i]+a.g[1][1][i]*b.g[1][1][i])%H;
  }

  for(int i=0;i<2;i++) for(int j=0;j<2;j++) NTT(c.g[i][j],P,-1),c.g[i][j].resize(K);
  return c;
}

但是,我们应该去研究一些倍增 NTT。

考虑能不能从 fafb 转移到 fa+b

发现是可以的

fa+b,i=j=0ifa,jfb,ji+j=0i1fa1,jfb1,ij1

很明显的卷积形式,表示出来就是

Fa+b(x)=Fa(x)Fb(x)+xFa1(x)Fb1(x)

于是我们尝试从 n2n,也就是

{F2n(x)=Fn(x)2+xFn1(x)2F2n1(x)=Fn(x)Fn1(x)+xFn1(x)Fn2(x)F2n2(x)=Fn1(x)2+Fn2(x)2

那么每次我们只需要维护 Fn(x),Fn1(x),Fn2(x) 就可以通过上述转移转移到 F2n(x),F2n1(x),F2n2(x)

这就是 倍增 NTT 了,直接做还是 O(klognlogk)。没有实现(应该和上述代码类似),咕。


另外,这里还有一种 O(klogk) 的高级做法。

观察上面的式子

Fi(x)=(x+1)Fi1(x)+xFi2(x)

发现是一个递推方程的形式,所以我们考虑用 OGF 去优化它。

也就是设

G(z)=F0+F1z+F2z2+

于是根据递推形式可以得到

G(z)=(x+1)G(z)z+xG(z)z2+1

转化出来可以得到

G(z)=11(1+x)zxz2

我们需要求出 Fn,于是想把这个式子展开,展开成

u1Az+v1Bz

的形式。


于是直接列方程

{A+B=1+xAB=x{A=(1+x)+(1+x)2+4x2B=(1+x)(1+x)2+4x2

带入 u,v 的方程

{u+v=1uB+vA=0{u=AABv=BAB

那么我们需要求的 [zn]G(z) 就可以根据 G(z) 的系数得到了

G(z)=un=0Anzn+vn=0Bnzn[zn]G(z)=uAn+vBn=An+1Bn+1AB

注意到 Bn+10(modxn+1),所以其实答案式子是

An+1AB

最后的式子就可以以 O(klogk) 的时间求得了。代码

这是一个 poly 版的板子!

Code
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define poly vector<ull>
#define pb push_back

const int N=2e6+5,H=998244353;
int n,m,U[N],inv[N],k;
ull w[N];

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int dec(int a,int b){return a<b?a-b+H:a-b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

const int g=3,ig=qpow(g);

void init(){
  inv[1]=1;
  for(int i=2;i<N;i++) inv[i]=mul(inv[H%i],dec(0,H/i));
}

void Init(int n){
  for(int i=0;i<n;i++) U[i]=(U[i>>1]>>1)|((i&1)?n>>1:0);
}

void Cro(poly &F,poly G){
  for(int i=0;i<(int)F.size();i++) F[i]=F[i]*G[i]%H;
}

void Clr(poly &F,int l,int r){
  for(int i=l;i<r;i++) F[i]=0;
}

poly Dlt(poly F){
  for(int i=1;i<(int)F.size();i++) F[i-1]=1ll*i*F[i]%H;
  F.pop_back();return F;
}

poly Sum(poly F){
  F.pb(0);
  for(int i=(int)F.size()-1;i>0;i--) F[i]=1ll*inv[i]*F[i-1]%H;
  F[0]=0;return F;
}

void NTT(poly &p,int len,int op){
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  ull u,tmp;w[0]=1;
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=qpow(op==1?g:ig,(H-1)/k);
	for(int i=1;i<=lst;i++) w[i]=w[i-1]*u%H;
	for(int i=0;i<len;i+=k)
	  for(int l=i;l<i+lst;l++)
	    tmp=p[l+lst]*w[l-i]%H,p[l+lst]=p[l]-tmp+H,p[l]+=tmp;
    if(k==(1<<16)) for(int i=0;i<len;i++) p[i]%=H;
  }
  for(int i=0;i<len;i++) p[i]%=H;
  if(op==-1){
    int I=qpow(len);
	for(int i=0;i<len;i++) p[i]=1ll*I*p[i]%H;
  }
}

poly operator *(poly F,poly G){
  int n=(int)F.size(),m=(int)G.size(),P;
  for(P=1;P<n+m;P<<=1);
  Init(P),F.resize(P),G.resize(P);

  NTT(F,P,1),NTT(G,P,1),Cro(F,G),NTT(F,P,-1);

  F.resize(n+m-1);return F;
}

poly Inv(poly F){
  int n=(int)F.size(),P;
  for(P=1;P<n;P<<=1);
  poly A,B,C;

  A={(ull)qpow(F[0])};
  for(int k=2;k<=P;k<<=1){
	B=F,B.resize(k),C=A,A.resize(k),Init(k);

	NTT(B,k,1),NTT(A,k,1),Cro(B,A),NTT(B,k,-1);
	Clr(B,0,k>>1),NTT(B,k,1),Cro(B,A),NTT(B,k,-1);

	for(int i=0;i<(k>>1);i++) A[i]=C[i];
	for(int i=(k>>1);i<k;i++) A[i]=(H-B[i])%H;
  }
  A.resize(n);
  return A;
}

poly Ln(poly F){
  poly A=Inv(F);int n=(int)F.size();
  F=Dlt(F),F=F*A;
  F=Sum(F),F.resize(n);
  return F;
}

poly Exp(poly F){
  int n=(int)F.size(),P;
  for(P=1;P<n;P<<=1);
  poly A,B;A={1},F.resize(P);

  for(int k=2;k<=P;k<<=1){
	B=A,B.resize(k),B=Ln(B);
	for(int i=0;i<k;i++) B[i]=(F[i]-B[i]+H)%H;
	B[0]=(B[0]+1)%H,A=A*B,A.resize(k);
  }
  A.resize(n);return A;
}

poly Qpow(poly F,int k){
  F=Ln(F);
  for(int i=0;i<(int)F.size();i++) F[i]=1ll*k*F[i]%H;
  return Exp(F);
}

poly Sqrt(poly F){
  int P,I=qpow(2),n=(int)F.size();
  for(P=1;P<n;P<<=1);
  poly A,B,C;

  A={1};

  for(int k=2;k<=P;k<<=1){
    A.resize(k),B=Inv(A);
	C=F,C.resize(k),B=B*C;
	for(int i=0;i<(k>>1);i++) B[i]=(B[i]+A[i])%H;
	for(int i=0;i<k;i++) A[i]=1ll*I*B[i]%H;
  }
  A.resize(n);return A;
}

poly A,C;

int main(){
  /*2024.7.20 H_W_Y CF755G PolandBall and Many Other Balls 多项式优化 dp*/
  init();int I=qpow(2);
  cin>>n>>m;k=m,m=min(m+1,n+1);
  C={1,6,1};C.resize(m);C=Sqrt(C);
  A=C,A[0]=(A[0]+1)%H,A[1]=(A[1]+1)%H;
  for(int i=0;i<m;i++) A[i]=1ll*I*A[i]%H;
  A=Qpow(A,n+1)*Inv(C);

  A.resize(m),A.resize(k+1);
  for(int i=1;i<=k;i++) cout<<A[i]<<' ';
  return 0;
}

P4173 残缺的字符串 - FFT 在字符串匹配中的应用

很久很久以前,在你刚刚学习字符串匹配的时候,有两个仅包含小写字母的字符串 AB,其中 A 串长度为 mB 串长度为 n。可当你现在再次碰到这两个串时,这两个串已经老化了,每个串都有不同程度的残缺。

你想对这两个串重新进行匹配,其中 A 为模板串,那么现在问题来了,请回答,对于 B 的每一个位置 i,从这个位置开始连续 m 个字符形成的子串是否可能与 A 串完全匹配?

1mn3×105

他都带通配符了,你就顺从它吧。


显然带着通配符,字符串的函数就不太能用了。

那怎么办呢?

考虑定义配对函数

Cx,y=[AxBy]2AxBy

当前字母是通配符时 Ax=0,容易发现,当且仅当 x,y 匹配时 C0

考虑设 Pk 表示 k 位置为结尾的匹配情况,于是 k 位置时合法的,当且仅当

Pk=i=0m1[AiBkm+i+1]2AiBkm+i+1=0

这并不好看。考虑按照之前的套路,我们令 Si=Am1i,也就可以得到

Pk=i=0m1[SiBki]2SiBki

这怎么这么像卷积形式?直接展开,可以得到

Pk=i=0kSi3Bki+i=0kSiBki32i=0kSi2Bk12

这三个部分都可以直接卷积处理了。

于是就做完了,FFT 即可。代码

Code
const int N=2.1e6+5;
const ld PI=acos(-1);
int n,m,U[N];
ll a[N],b[N],c[N],d[N],e[N];
cp F[N];
string s;
vector<int> ans;

void FFT(cp *p,int len,int op){}
void Mul(ll *a,ll *b,int n,int m){}

int main(){
  /*2024.7.11 H_W_Y P4173 残缺的字符串 FFT*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>m>>n>>s;
  for(int i=0;i<m;i++){
	if(s[i]=='*') a[m-i-1]=0;
	else a[m-i-1]=s[i]-'a'+1;
  }
  for(int i=0;i<m;i++) b[i]=a[i]*a[i],c[i]=b[i]*a[i];
  cin>>s;
  for(int i=0;i<n;i++){
	if(s[i]=='*') d[i]=0;
	else d[i]=s[i]-'a'+1;
	e[i]=d[i];
  }
  Mul(c,d,m,n);
  for(int i=0;i<n;i++) d[i]*=e[i];
  Mul(b,d,m,n);
  for(int i=0;i<n;i++) d[i]*=e[i];
  Mul(a,d,m,n);

  for(int i=m-1;i<n;i++)
    if(a[i]+c[i]-2*b[i]==0) ans.pb(i-m+2);
  
  cout<<(int)ans.size()<<'\n';
  for(auto i:ans) cout<<i<<' ';
  return 0;
}

CF528D Fuzzy Search - FFT 在字符串匹配中的应用

给出一个门限值 k 和两个只包含 AGCT 四种字符的基因串 ST。现在你要找出在下列规则中 TS 中出现了几次。

TS 的第 i 个位置中出现,当且仅当把 T 的首字符和 S 的第 i 个字符对齐后,T 中的每一个字符能够在 S 中找到一个位置偏差不超过 k 的相同字符。

即对于所有的 j[1,|T|],都存在一个 p[1,|S|] 使得 |(i+j1)p|kSp=Tj

例如 k=1 时,ACAT 出现在 AGCAATTCAT2 号, 3 号和 6 号位置。 (编号从 1 开始。)

1nm2105 , 0k2105

居然切了。


考虑按照上一题的思路,我们希望找到一个匹配函数来表示是否匹配这件事情。

由于 DNA 中只有四种碱基,所以我们考虑把它们每一个单独拿出来,构造一下两个序列

  • A:若 Ai=1 则表示 i 位置合法,即存在一个和 i 距离不超过 k 的碱基;反之则 Ai=0
  • B:若 Bi=1 则表示 Tm1i 是当前的碱基,反之则不是。

于是 k 位置为结尾的字符串匹配当且仅当

Pk=i=0kBiAki=0

就可以用 FFT 直接完成了。代码

Code
const int N=2.1e6+5;
const ld PI=acos(-1);
int n,m,U[N],k,c[N],ans=0;
string s,t;
cp F[N];
ll a[N],b[N];

void FFT(cp *p,int len,int op){}
void Mul(ll *a,ll *b,int n,int m){}

void SOLVE(char ch){
  for(int i=0;i<m;i++) b[m-i-1]=t[i]==ch;
  for(int i=0;i<n;i++) a[i]=k+1;

  int lst=-1;
  for(int i=0;i<n;i++){
	if(s[i]==ch) lst=i;
	if(lst>=0) a[i]=min(a[i],(ll)i-lst);
  }
  lst=n+1;
  for(int i=n-1;i>=0;i--){
	if(s[i]==ch) lst=i;
	if(lst<n) a[i]=min(a[i],(ll)lst-i);
  }
  for(int i=0;i<n;i++) a[i]=(a[i]>k);

  Mul(a,b,n,m);

  for(int i=m-1;i<n;i++) if(!a[i]) c[i-m+1]++;
}

int main(){
  /*2024.7.11 H_W_Y CF528D Fuzzy Search FFT*/
  cin>>n>>m>>k>>s>>t;
  SOLVE('A'),SOLVE('T'),SOLVE('G'),SOLVE('C');
  for(int i=0;i<n;i++) if(c[i]==4) ++ans;
  cout<<ans;
  return 0;
}

CF622F The Sum of the k-th Powers - 拉插

(i=1nik)mod(109+7)

1n109,0k106

有些人拉插喜欢把 xi 当成 yi 算。


首先,为什么这是一个 k+1 次多项式?

一个序列做 p 阶差分变成常数,当且仅当这是一个 p 次多项式。

因为每一次差分会让次数减少 1,这个结论也是可以用牛顿级数证明的。

由于我们把这个式子差分之后是 1k,2k,,nkk 次多项式,所以原式是一个 k+1 次多项式。


于是似乎就做完了,直接拉插

i=1k+2yijinjij

这东西稍微预处理一下就容易完成了。代码

Code
#include <bits/stdc++.h>
using namespace std;

const int N=1e6+5,H=1e9+7;
int n,m,suf[N],ifac[N],ans=0;

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int dec(int a,int b){return a<b?a-b+H:a-b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

void init(){
  ifac[m+1]=suf[m+3]=1;
  for(int i=1;i<=m+1;i++) ifac[m+1]=mul(ifac[m+1],i);
  ifac[m+1]=qpow(ifac[m+1]);
  for(int i=m+1;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
  for(int i=m+2;i>=1;i--) suf[i]=mul(suf[i+1],dec(n,i));
}

int main(){
  /*2024.7.11 H_W_Y CF622F The Sum of the k-th Powers 插值*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>m;
  init();
  for(int i=1,nw=0,a=1,tmp;i<=m+2;i++){
    add(nw,qpow(i,m));
    tmp=mul(nw,mul(mul(a,suf[i+1]),mul(ifac[i-1],ifac[m+2-i])));
    if((m+2-i)&1) del(ans,tmp);
	else add(ans,tmp);
	a=mul(a,dec(n,i));
  }
  cout<<ans<<'\n';
  return 0;
}

CF623E Transforming Sequence - 倍增多模 NTT

对于一个整数序列 {a1,a2,,an},我们定义它的变换为 {b1,b2,,bn},其中 bi=a1|a2||ai,其中 | 表示二进制按位或运算。

现在求有多少个长为 n 的序列 a,满足 1ai2k1,使得它的变换 b严格单调递增的,对 109+7 取模。

1n10181k3×104

发现多模 NTT 写成 vector 形式其实挺好写的。


首先,容易发现 n>k 的答案是 0,于是接下来我们默认 nk 了。

写出 dp 方程,设 fi,j 表示前 i 个数有 j 位为 1 的方案数,我们不关心是哪 j 位,也就是说与 k 是无关的。

那么可以得到

fi,j=p=1jfi1,jp2jp(jp)

看起来就非常卷积。

于是考虑倍增 NTT,考虑知道 fa,fb 时如何推出 fa+b,容易列出式子

fa+b,j=p=1jfa,pfb,jp(jp)2bpfa+b,jj!=p=0jfa,p2bpp!fb,jp(jp)!

就可以直接用多模 NTT 转移了。/qiang

时间复杂度 O(klogklogn)代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ull unsigned long long
#define bi __int128
#define poly vector<ull>

const int N=2e5+5,H=1e9+7;
int K,U[N],P,ifac[N],fac[N],ans=0;
ull w[N],n;
map<int,poly> mp;

int qpow(int a,int b,int p){
  int res=1;
  while(b){if(b&1) res=1ll*res*a%p;a=1ll*a*a%p,b>>=1;}
  return res;
}

void NTT(poly &a,int op,int p){
  int len=P;
  for(int i=0;i<len;i++) if(i<U[i]) swap(a[i],a[U[i]]);
  ull u,tmp,ig=qpow(3,p-2,p);w[0]=1;
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=qpow(op==1?3:ig,(p-1)/k,p);
	for(int i=1;i<=lst;i++) w[i]=w[i-1]*u%p;
	for(int i=0;i<len;i+=k)
	  for(int l=i;l<i+lst;l++)
	    tmp=w[l-i]*a[l+lst]%p,a[l+lst]=a[l]-tmp+p,a[l]+=tmp;
	if(k==(1<<16)) for(int i=0;i<len;i++) a[i]%=p;
  }
  for(int i=0;i<len;i++) a[i]%=p;
  if(op==-1){
	int I=qpow(len,p-2,p);
	for(int i=0;i<len;i++) a[i]=a[i]*I%p;
  }
}

poly Mul(poly a,poly b,int p){
  a.resize(P),b.resize(P);
  NTT(a,1,p),NTT(b,1,p);
  for(int i=0;i<P;i++) a[i]=a[i]*b[i]%p;
  NTT(a,-1,p);
  return a;
}

poly operator *(poly a,poly b){
  poly res[3];
  a.resize(P),b.resize(P);
  for(int i=P>>1;i<P;i++) a[i]=b[i]=0;
  bi p[3]={998244353,469762049,167772161},M=p[0]*p[1]*p[2];
  for(int i=0;i<3;i++) res[i]=Mul(a,b,p[i]);
  for(int i=0;i<P;i++){
	bi v=0;
	for(int j=0;j<3;j++) v+=res[j][i]*(M/p[j])*qpow((M/p[j])%p[j],p[j]-2,p[j]);
	a[i]=v%M%H;
  }
  return a;
}

poly solve(int n){
  poly h(P,0);
  if(n==1){
    for(int i=1;i<=K;i++) h[i]=1;
	return mp[n]=h;
  }
  if(mp.count(n)) return mp[n];
  int a=n>>1,b=n-a,pw=qpow(2,b,H);
  poly f=solve(a),g=solve(b);

  for(int i=1,nw=1;i<=K;i++) f[i]=1ll*f[i]*ifac[i]%H*(nw=1ll*nw*pw%H)%H;
  for(int i=1;i<=K;i++) g[i]=1ll*g[i]*ifac[i]%H;

  h=f*g;
  for(int i=1;i<=K;i++) h[i]=1ll*h[i]*fac[i]%H;

  return mp[n]=h;
}

void init(){
  for(P=1;P<=(K<<1);P<<=1);
  for(int i=0;i<P;i++) U[i]=(U[i>>1]>>1)|((i&1)?P>>1:0);

  fac[0]=1;
  for(int i=1;i<=P;i++) fac[i]=1ll*fac[i-1]*i%H;
  ifac[P]=qpow(fac[P],H-2,H);
  for(int i=P;i>=1;i--) ifac[i-1]=1ll*ifac[i]*i%H;
}

int binom(int n,int m){
  if(m<0||n<m) return 0;
  return 1ll*fac[n]*ifac[m]%H*ifac[n-m]%H;
}

int main(){
  /*2024.7.11 H_W_Y CF623E Transforming Sequence 倍增多模 NTT*/
  ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
  cin>>n>>K;init();
  if(n>(ull)K) cout<<"0\n",exit(0);
  
  poly res=solve(n);

  for(int i=n;i<=K;i++) (ans+=1ll*binom(K,i)*res[i]%H)%=H;
  cout<<ans<<'\n';

  return 0;
}

P3723 [AH2017/HNOI2017] 礼物 - FFT

我的室友最近喜欢上了一个可爱的小女生。马上就要到她的生日了,他决定买一对情侣手环,一个留给自己,一个送给她。每个手环上各有 n 个装饰物,并且每个装饰物都有一定的亮度。

但是在她生日的前一天,我的室友突然发现他好像拿错了一个手环,而且已经没时间去更换它了!他只能使用一种特殊的方法,将其中一个手环中所有装饰物的亮度增加一个相同的非负整数 c。并且由于这个手环是一个圆,可以以任意的角度旋转它,但是由于上面装饰物的方向是固定的,所以手环不能翻转。需要在经过亮度改造和旋转之后,使得两个手环的差异值最小。

在将两个手环旋转且装饰物对齐了之后,从对齐的某个位置开始逆时针方向对装饰物编号 1n,其中 n 为每个手环的装饰物个数, 第 1 个手环的 i 号位置装饰物亮度为 xi,第 2 个手环的 i 号位置装饰物亮度为 yi,两个手环之间的差异值为(参见输入输出样例和样例解释):

i=1n(xiyi)2

麻烦你帮他计算一下,进行调整(亮度改造和旋转),使得两个手环之间的差异值最小,这个最小值是多少呢?

1n50000, 1aim100

居然自己做出来了 /jk


首先考虑把 b 倍长一次,然后考虑移动 i 位,那么就有

Pi=j=0n1(ajbj+i)2

根据之前字符串匹配的经验,和翻转的套路,这东西一看就非常卷积。

于是我们设

Pi=P2n1i,bi=b2n1i

那么

Pi=j=0n1aj2+j=0n1bij22j=0iajbij

简直就是卷积。于是可以 FFT 直接做。

而假设变化量是 c,我们直接在 b 的每一项加上一个 c,拆下括号发现是可以 O(1) 计算的。

于是就做完了,只需要一次 FFT 即可。代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ld long double
#define ll long long
#define cp complex<ld>

const int N=4.2e6+5;
const ld PI=acos(-1);
int n,m,U[N];
cp F[N];
ll a[N],b[N],suma,sumb,sa,sb,ans=(ll)1e18;

void FFT(cp *p,int len,int op){}
void Mul(ll *a,ll *b,int n,int m){}

int main(){
  /*2024.7.11 H_W_Y P3723 [AH2017/HNOI2017] 礼物 FFT*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>m;
  for(int i=0;i<n;i++) cin>>a[i];
  for(int i=0;i<n;i++) cin>>b[n-i-1],b[2*n-i-1]=b[n-i-1];
  for(int i=0;i<n;i++) sa+=a[i],sb+=b[i],suma+=a[i]*a[i],sumb+=b[i]*b[i];

  Mul(a,b,n,(n<<1));

  for(int i=0;i<(n<<1);i++) a[i]=suma+sumb-2ll*a[i];

  for(ll c=-m;c<=m;c++)
    for(int i=n;i<(n<<1);i++)
	  ans=min(ans,a[i]+1ll*n*c*c+2ll*c*(sb-sa));
  cout<<ans;
  return 0;
}

P4841 [集训队作业2013] 城市规划 - Ln

你需要求出 n 个点的简单 (无重边无自环) 有标号无向连通图数目。

由于这个数字可能非常大, 你只需要输出方案数对 1004535809 ( 479×221+1 ) 取模即可。

n130000

有些人板子写挂调了一万年。


考虑经典结论,若 A 中包含若干个 B,则 A^(x)=exp(B^(x))

其中 A^(x) 表示 A 的指数生成函数,B^(x) 同理。

于是设 A^(x)=n=02(n2)n!xn,即带标号的图个数,于是答案 B^(x)=ln(A^(x)),结束!

时间复杂度 O(nlogn)代码

Code
int main(){
  /*2024.7.16 H_W_Y P4841 [集训队作业2013] 城市规划 Ln*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n;init();++n;
  for(int i=0;i<n;i++) F[i]=1ll*qpow(2,(1ll*i*(i-1)/2))*ifac[i]%H;
  Ln(F,n);
  cout<<1ll*F[n-1]*fac[n-1]%H<<'\n';
  return 0;
}

P4389 付公主的背包 - Exp

这个背包最多可以装 105 大小的东西

付公主有 n 种商品,她要准备出摊了

每种商品体积为 vi,都有无限件

给定 m,对于 s[1,m],请你回答用这些商品恰好装 s 体积的方案数

1n,m1051vim

1.17s 略胜一筹。


考虑设 ai 表示体积为 i 的物品种类数,于是我们有生成函数

1im(1+xi+x2i+)ia=1im(11xi)ai

不知道怎么做?经典操作:先去 Ln,再去 Exp,可以得到

1im(11xi)ai=exp(i=1mailn(1xi))=exp(i=1maij=1mixijj)

后面的式子是调和级数形式,只有 O(nlogn) 项,可以暴力求和。

于是最后再做一次 exp 就可以了,时间复杂度 O(mlogm)。代码。

Code
int main(){
  /*2024.7.16 H_W_Y P4389 付公主的背包 多项式 Exp*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>m;++m;init();
  for(int i=1,x;i<=n;i++) cin>>x,a[x]++;
  for(int i=1;i<=m;i++) if(a[i]){
	for(int j=1;j*i<=m;j++)
	  (F[i*j]+=mul(inv[j],a[i]))%=H;
  }
  Exp(F,m);
  for(int i=1;i<m;i++) cout<<F[i]<<'\n';
  return 0;
}

UOJ182 a^-1 + b problem - 多点求值

首先炮身的屏幕上显示了 n 个数,接着在模 998244353(7×17×223+1) 一个质数意义下,给出了 m 条指令,每一条指令都是下列两种之一:

  1. 所有数加上某一个数。
  2. 所有数都变成原来的逆元。(保证这时所有数的逆元都存在)

在每个指令完成之后,你要回答当前屏幕上所有数的和。

跳蚤国王思索了片刻,发现这个问题奥妙重重,于是他让你——这附近最著名的民间科学家来帮他解决这个难题。

1n105,1m60000

多项式写成 poly 还是要优秀一点/kk


我们设初始的序列为 x1,x2,,xn,那么第 i 次操作后,xj 一定形如

aixj+bicixj+di

而这里的 ai,bi,ci,di 是可以线性求出来的。

进而,我们可以将这个式子化简,根据初中数学的套路,该式可以被化简成

1xj+ti×ki+ei

我们发现求和之后相当于是去求

F(t)=i=1n1xi+t

t=ti 时的值,那么这就是一个多点求值了。


而这个式子不是标准的多项式,看到它应该想到 求导

对式子进行通分,那么分母就是

i=1n(xi+t)

然而分子——是分母的 导数!这是可以用链式法则证明的。

于是就做完了,直接多点求值即可。代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define poly vector<ull>
#define pb push_back

const int N=5e5+5,H=998244353;
int n,m,U[N],x[N],k[N],a,b,c,d,e[N],sum,as[N];
ull w[N];
poly F,G,t;

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int dec(int a,int b){return a<b?a-b+H:a-b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

ll qpow(ll a,ll b=H-2){
  ll res=1;
  while(b){if(b&1) res=res*a%H;a=a*a%H,b>>=1;}
  return res;
}

const int g=3,ig=qpow(g);

void Init(int n){
  for(int i=0;i<n;i++) U[i]=(U[i>>1]>>1)|((i&1)?n>>1:0);
}

void Cro(poly &F,poly G,int n){
  for(int i=0;i<n;i++) F[i]=F[i]*G[i]%H;
}

void Clr(poly &F,int l,int r){
  for(int i=l;i<r;i++) F[i]=0;
}

void Rev(poly &F,int l,int r){
  reverse(F.begin()+l,F.begin()+r);
}

poly Dlt(poly F){
  for(int i=1;i<(int)F.size();i++) F[i-1]=1ll*i*F[i]%H;
  F.pop_back();return F;
}

void NTT(poly &p,int len,int op){
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  ull u,tmp;w[0]=1;
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=qpow(op==1?g:ig,(H-1)/k);
	for(int i=1;i<=lst;i++) w[i]=w[i-1]*u%H;
	for(int i=0;i<len;i+=k)
	  for(int l=i;l<i+lst;l++)
	    tmp=w[l-i]*p[l+lst]%H,p[l+lst]=p[l]-tmp+H,p[l]+=tmp;
	if(k==(1<<16)) for(int i=0;i<len;i++) p[i]%=H;
  }
  for(int i=0;i<len;i++) p[i]%=H;
  if(op==-1){
	int I=qpow(len);
	for(int i=0;i<len;i++) p[i]=1ll*I*p[i]%H;
  }
}

poly operator *(poly F,poly G){
  int n=(int)F.size(),m=(int)G.size();
  int P;for(P=1;P<n+m;P<<=1);
  Init(P),F.resize(P),G.resize(P);
  NTT(F,P,1),NTT(G,P,1),Cro(F,G,P),NTT(F,P,-1);
  F.resize(n+m-1);return F;
}

void Inv(poly &F){
  int n=(int)F.size();
  int P;for(P=1;P<n;P<<=1);
  poly A,B,C;

  A.pb(qpow(F[0]));
  for(int k=2;k<=P;k<<=1){
	B=F,B.resize(k),C=A,A.resize(k),Init(k);

	NTT(B,k,1),NTT(A,k,1),Cro(B,A,k),NTT(B,k,-1);
	Clr(B,0,k>>1),NTT(B,k,1),Cro(B,A,k),NTT(B,k,-1);

	for(int i=0;i<(k>>1);i++) A[i]=C[i];
	for(int i=(k>>1);i<k;i++) A[i]=(H-B[i])%H;
  }
  A.resize(n);F=A;
}

poly operator % (poly F,poly G){
  int n=(int)F.size(),m=(int)G.size();
  if(n<m) return F;
  poly A,B;
  int P=n-m+1;
  Rev(F,0,n),A=F,A.resize(P),Rev(F,0,n);
  Rev(G,0,m),B=G,B.resize(P),Rev(G,0,m);

  Inv(B),B=A*B,Rev(B,0,P),B.resize(P),B=B*G;
  for(int i=0;i<m-1;i++) G[i]=(F[i]-B[i]+H)%H;
  G.resize(m-1);return G;
}

#define mid ((l+r)>>1)
#define lc p<<1
#define rc p<<1|1
#define lson l,mid,lc
#define rson mid+1,r,rc

poly build0(int l,int r,int p){
  if(l==r) return {(ull)x[l],1};
  return build0(lson)*build0(rson);
}

poly Q[N<<2];

void build(int l,int r,int p,poly &a){
  if(l==r) return Q[p]={(H-a[l])%H,1},void();
  build(lson,a),build(rson,a),Q[p]=Q[lc]*Q[rc];
}

void solve(int l,int r,int p,poly F,poly G){
  if(l==r){
	if(e[l]==-1) return;
	as[l]=mul(G[0],qpow(F[0]));
	return;
  }
  solve(lson,F%Q[lc],G%Q[lc]),solve(rson,F%Q[rc],G%Q[rc]);
}

#define H_W_Y

int main(){
  #ifdef H_W_Y
  freopen("1.in","r",stdin);
  freopen("1.out","w",stdout);
  #endif

  /*2024.7.16 H_W_Y UOJ #182. 【UR #12】a^-1 + b problem 多项式多点求值*/
  ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);
  cin>>n>>m;t.resize(m);
  for(int i=0;i<n;i++) cin>>x[i],add(sum,x[i]);
  a=d=1,b=c=0;

  for(int i=0,op,v;i<m;i++){
	cin>>op;
	if(op==1) cin>>v,add(a,mul(c,v)),add(b,mul(d,v));
    else swap(a,c),swap(b,d);

	if(c==0){
      e[i]=-1,t[i]=0;
	  k[i]=adc(mul(sum,mul(a,qpow(d))),mul(mul(n,b),qpow(d)));
	}else{
      e[i]=mul(n,mul(a,qpow(c)));
	  t[i]=mul(d,qpow(c));
	  k[i]=mul(dec(mul(b,c),mul(a,d)),qpow(mul(c,c)));
	}
  }

  F=build0(0,n-1,1),G=Dlt(F);
  build(0,m-1,1,t),solve(0,m-1,1,F%Q[1],G%Q[1]);

  for(int i=0;i<m;i++){
	if(e[i]==-1) cout<<k[i]<<'\n';
	else cout<<adc(e[i],mul(k[i],as[i]))<<'\n';
  }

  return 0;
}

P5450 [THUPC2018] 淘米神的树 - 多点求值

可爱的 Tommy 有一棵树。这棵树上只有点 ab 是黑色,其它的点都是白色。

每次,Tommy 可以将一个黑色 p 的点染成红色,然后把和 p 相邻的所有白色的点染成黑色。最后,所有的点都会被染成红色。

设第 i 个点是第 ti 个被染成红色的,那么 ti 是一个 1n 的排列。Tommy 希望你帮他求出,有多少种不同的 ti

a,bn234,567

一发过编过样例过题(怎么根 pkusc 场上一样),然后喜提最劣解(毕竟是取模)。


考虑不在 ab 路径上的子树,那么每一次一定是最先染根,这就是一个满足树形拓扑序的排列问题。

和 6.29 T3 一样,结论:随机排列满⾜树形拓扑序的概率是⼦树⼤小的倒数乘积

证明可以去看那场联考的总结,写在里面了(7月训练记录)。


于是按照这个思路想下去,我们关心地其实就是 ab 路径上点染色的顺序。

由于有两个点,我们尝试把它合成一个,也就是在建立一个虚点 s,向 a,b 分别连无向边。

那么现在就变成了基环树,我们只关心环上的点的顺序,而对于不在环上的点,它子树的贡献是一定的。

先给环上点标号,s0,其余依次从 1k,其中 k 表示环上的点的个数。


枚举最后一个被染上的点?

怎么搞都是 n2 的 dp,怎么办呢?

由于是基环树——我们能不能考虑 断掉 一条边让他变成树?

发现是可以的,我们枚举基环树上的一条边 rr+1,断掉它表示最后一个点要么是 r 要么是 r+1,于是最后需要乘上 12

我们再钦定每个点的权值是 ai 表示环上第 i 个点,其非环上的子树大小之和加一;Siai 的前缀和。

那么答案就是

12C(i=0k1j=1i(SiSj1)j=i+1k(SjSi))

其中 C 是非环上子树的贡献。

于是现在我们研究的其实是

j=1i(SiSj1)j=i+1k(SjSi)

想到 导数。发现这个式子其实是函数

Fi(x)=j=0k(xSj)xSi

xSi 时的值再乘上一个 (1)ki

根据 洛必达法则,我们对上下同时求导,也就是求

F(x)=(j=0k(xSj))

x=Si 处的值。


由于我们要同时求每个 x=Si 的值,所以这是需要用多项式多点求值完成的。

于是就做完了,最终的答案式子就是

ans=12C(i=0k1F(Si)(1)ki)

时间复杂度 O(nlog2n)代码

Code
/*上面的代码和上一道题几乎一样*/
poly Q[N<<2];

void build(int l,int r,int p){
  if(l==r) return Q[p]={(ull)(H-s[l])%H,1},void();
  build(lson),build(rson),Q[p]=Q[lc]*Q[rc];
}

void solve(int l,int r,int p,poly F){
  if(l==r) return as[l]=F[0],void();
  solve(lson,F%Q[lc]),solve(rson,F%Q[rc]);
}

void dfs(int u){
  sz[u]=1,f[u][0]=fa[u],dep[u]=dep[fa[u]]+1;
  for(int i=1;f[f[u][i-1]][i-1];i++) f[u][i]=f[f[u][i-1]][i-1];
  for(auto v:G[u]) if(v!=fa[u]) fa[v]=u,dfs(v),sz[u]+=sz[v];
}

int LCA(int u,int v){
  if(dep[u]<dep[v]) swap(u,v);
  for(int i=18;i>=0;i--) if(dep[f[u][i]]>=dep[v]) u=f[u][i];
  if(u==v) return u;
  for(int i=18;i>=0;i--) if(f[u][i]!=f[v][i]) u=f[u][i],v=f[v][i];
  return f[u][0];
}

void dfs2(int u,int pre){
  sz[u]=1;
  for(auto v:G[u]) if(!vis[v]&&v!=pre) dfs2(v,u),sz[u]+=sz[v];
  ans=mul(ans,qpow(sz[u]));
}

void SOLVE(){
  ans=0,dfs(1);
  int lca=LCA(a,b),u=a;
  while(u!=lca){
    p[++ct]=u,vis[u]=true;
	for(auto v:G[u])
	  if(v!=fa[u]&&v!=p[ct-1])
	    s[ct]+=sz[v];
	++s[ct];
	u=fa[u];
  }
  int nw=ct+dep[b]-dep[lca]+1;u=b;
  while(u!=lca){
	p[nw]=u,vis[u]=true;
	for(auto v:G[u])
	  if(v!=fa[u]&&v!=p[nw+1])
	    s[nw]+=sz[v];
	++s[nw];
	u=fa[u],--nw;
  }
  ++ct;
  p[ct]=lca,vis[lca]=true;
  for(auto v:G[lca]) if(v!=fa[lca]&&v!=p[ct-1]&&v!=p[ct+1]) s[ct]+=sz[v];
  s[ct]+=n-sz[lca]+1;
  ct+=dep[b]-dep[lca];

  for(int i=1;i<=ct;i++) s[i]+=s[i-1];

  build(0,ct,1),F=Dlt(Q[1]);
  solve(0,ct,1,F);

  for(int i=0;i<=ct;i++){
	if((ct-i)&1) del(ans,qpow(as[i]));
    else add(ans,qpow(as[i]));
  }
  for(int i=1;i<=ct;i++) for(auto v:G[p[i]]) if(!vis[v]) dfs2(v,p[i]);

  ans=mul(ans,qpow(2));
  for(int i=1;i<=n;i++) ans=mul(ans,i);
}

int main(){
  /*2024.7.17 H_W_Y P5450 [THUPC2018] 淘米神的树 多项式多点求值*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  
  cin>>n>>a>>b;
  for(int i=1,u,v;i<n;i++) cin>>u>>v,G[u].pb(v),G[v].pb(u);

  SOLVE();

  cout<<ans<<'\n';
  return 0;
}

*降幂引理

根据名字就可以知道,这个引理用来解决指数 %k 的问题。

也就是我们现在有一个多项式 A(x),令

B(x)=i=0naiximodk

n 次多项式变成 k 次多项式。

然而这两个多项式满足

A(ωk0)A(ωkk1)=B(ωk0)B(ωkk1)

即在 k 次单位根的位置的点值相同。


这是容易证明的,考虑直接把式子化开

A(ωkj)=aiωij=ai(wi%k)j=ai(ωj)i%k=B(ωkj)

证毕。


降幂引理单位根反演 本质 相同!!!在之后的例题中我们会发现这一点。


*Bluestein

也叫 Chirp Z-Transform。

任意模长的 DFT,不要求 n2k。/qiang

曾经我们将卷积 DFT/IDFT,现在我们尝试把 DFT/IDFT 转化成卷积。

给定多项式 A(x)=a0++an1xn1,求其在 ω0,,ωn1 的值,不要求 ω 是单位根(可以是任意数,只要这是个等比数列就行)。

b0bn1,其中 bm=k=0n1ak(ωm)k

多项式多点求值?太慢啦!


考虑重要的 trick

mk=12[m2+(k+1)2+(mk)2]

展开即可验证。

不能展开成 12[mk+k2(mk)2],在后面的推到中不保证 ω 有二次剩余,所以非常局限!!!

于是我们把这个式子带入之前的式子中

bm=ω12m2i=0n1akω12(k+1)212(mk)2

这不就是卷积形式?但是又不标准。


于是我们考虑一种非常简单的转化,令

Ai=aiω12(k+1)2,Bm+n=bmω12m2,Cj=ω12(jn)2

我们就有

Bm+n=i=0m+nAiCm+ni

这就是一个标准卷积形式了,直接 FFT 或者 NTT 即可。常数可能有点大。


对其进行一个推广。

A(x)kk 次单位根的值。

直接用降幂引理可以 A 转化成 k 次多项式,于是就可以用 Bluestein 做了。


P6800 【模板】Chirp Z-Transform

板板。

给定一个 n 项多项式 P(x) 以及 c,m,请计算 P(c0),P(c1),,P(cm1)。所有答案都对 998244353 取模。

1n,m106,0c,ai<998244353.


没有预处理,直接暴力做!代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ull unsigned long long
#define ll long long

const int N=4.2e6+5,H=998244353;
int n,m,c,U[N],a[N],k;
ull A[N],C[N],w[N];

ll qpow(ll a,ll b=H-2){
  ll res=1;a%=H,b%=(H-1),b=(b+H-1)%(H-1);
  while(b){if(b&1) res=1ll*res*a%H;a=1ll*a*a%H,b>>=1;}
  return res;
}

const int g=3,ig=qpow(g);

void NTT(ull *p,int len,int op){
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  ull u,tmp;w[0]=1;
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=qpow(op==1?g:ig,(H-1)/k);
	for(int i=1;i<=lst;i++) w[i]=w[i-1]*u%H;
	for(int i=0;i<len;i+=k)
	  for(int l=i;l<i+lst;l++)
	    tmp=w[l-i]*p[l+lst]%H,p[l+lst]=p[l]-tmp+H,p[l]+=tmp;
	if(k==(1<<16)) for(int i=0;i<len;i++) p[i]%=H;
  }
  for(int i=0;i<len;i++) p[i]%=H;
  if(op==-1){
	ull I=qpow(len);
	for(int i=0;i<len;i++) p[i]=p[i]*I%H;
  }
}

void Init(int n){
  for(int i=0;i<n;i++) U[i]=(U[i>>1]>>1)|((i&1)?n>>1:0);
}

void Cro(ull *F,ull *G,int n){
  for(int i=0;i<n;i++) F[i]=F[i]*G[i]%H;
}

void Mul(ull *F,ull *G,int n,int m){
  int P=1;
  for(P=1;P<n+m;P<<=1);
  Init(P),NTT(F,P,1),NTT(G,P,1),Cro(F,G,P),NTT(F,P,-1);
}

int main(){
  /*2024.7.15 H_W_Y P6800 【模板】Chirp Z-Transform Bluestein*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>c>>m;
  for(int i=0;i<n;i++) cin>>A[i],A[i]=1ull*A[i]*qpow(c,1ll*i*(i+1)/2ll)%H;
  n=max(n,m);
  for(int i=0;i<(n<<1);i++) C[i]=qpow(c,-(1ll*(i-n)*(i-n-1)/2ll));
  Mul(A,C,n,(n<<1));
  for(int i=n;i<n+m;i++) cout<<A[i]*qpow(c,1ll*(i-n)*(i-n-1)/2ll)%H<<' ';
  return 0;
}

P6828 任意模数 Chirp Z-Transform

套上 MTT 即可。

注意预处理单位根。代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ld long double
#define cp complex<ld>
#define ll long long
#define il inline

const int N=2.1e6+5,H=1e9+7,Bl=1<<15;
const ld PI=acos(-1);
int n,m,c,U[N],a[N],b[N],fac[N],ifac[N];
cp A[N],B[N],C[N],w1[N],w2[N];

il ll qpow(ll a,ll b=H-2){
  ll res=1;a%=H,b=(b%(H-1)+H-1)%(H-1);
  while(b){if(b&1) res=1ll*res*a%H;a=1ll*a*a%H,b>>=1;}
  return res;
}

il void FFT(cp *p,int len,int op){
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  cp nw,tmp,u,*w;
  for(int k=2,lst=1,del=len;k<=len;lst=k,k<<=1,del>>=1){
	for(int i=0;i<len;i+=k){
      w=(op==1?&w1[0]:&w2[0]);
	  for(int l=i;l<i+lst;l++)
	    tmp=*w*p[l+lst],p[l+lst]=p[l]-tmp,p[l]+=tmp,w+=del;
	}
  }
}

il void Init(int n){
  for(int i=0;i<n;i++) U[i]=(U[i>>1]>>1)|((i&1)?n>>1:0);
  for(int i=0;i<n;i++){
	w1[i]=cp(cos(i*PI/n),sin(i*PI/n));
	w2[i]=cp(cos(i*PI/n),-sin(i*PI/n));
  }
}

il void Mul(int *a,int *b,int n,int m){
  int P=1;
  for(P=1;P<n+m;P<<=1);
  Init(P);

  for(int i=0;i<n;i++) A[i]+=cp(a[i]/Bl,a[i]%Bl),C[i]+=cp(a[i]/Bl,-a[i]%Bl);
  for(int i=0;i<m;i++) B[i]+=cp(b[i]/Bl,b[i]%Bl);

  FFT(A,P,1),FFT(B,P,1),FFT(C,P,1);
  for(int i=0;i<P;i++) A[i]=A[i]*B[i],C[i]=C[i]*B[i];
  FFT(A,P,-1),FFT(C,P,-1);

  for(int i=0;i<P;i++){
	ll x=(ll)((A[i].real()+C[i].real())/P/2+0.5),y=(ll)(A[i].imag()/P+0.5),z=(ll)((C[i].real()-A[i].real())/P/2+0.5);
	a[i]=1ll*(1ll*x%H*Bl%H*Bl%H+1ll*y%H*Bl%H+z%H)%H;
  }
}

int main(){
  /*2024.7.15 H_W_Y P6800 【模板】Chirp Z-Transform Bluestein*/
//   freopen("1.err","r",stdin);
//   freopen("1.out","w",stdout);
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>c>>m;
  for(int i=0;i<n;i++) cin>>a[i];
  n=max(n,m);
  
  fac[0]=1;
  for(int i=1;i<=(n<<1);i++) fac[i]=1ll*fac[i-1]*c%H;
  ifac[n<<1]=qpow(fac[n<<1]);
  for(int i=(n<<1);i>=1;i--) ifac[i-1]=1ll*ifac[i]*c%H;

  for(int i=0,pro=1;i<n;i++) pro=1ll*pro*fac[i]%H,a[i]=1ll*a[i]*pro%H;

  for(int i=n,pro=1;i<(n<<1);i++) b[i]=pro,pro=1ll*pro*ifac[i-n]%H;
  for(int i=n-1,pro=ifac[1];i>=0;i--) b[i]=pro,pro=1ll*pro*ifac[n+1-i]%H;
  Mul(a,b,n,(n<<1));
  for(int i=n,pro=1;i<n+m;i++) cout<<1ll*a[i]*pro%H<<' ',pro=1ll*pro*fac[i-n]%H;
  return 0;
}

似乎有更好的拆分

mk=(m+k2)(m2)(k2)

这可以避免一些神秘的负数问题。我也不知道是什么问题。


*DFT 多路分治算法

不如 Bluestein 应用广,但是不需要 FFT。


我们先来考虑 n=3k 时如何做 DFT?

给定 f(x)=a0+a1x++an1xn1,求它在 ω0,,ωn1 的值。(ω=ωn,是单位根)


发现我们可以类似 2k 的做法,把整个式子分成三个部分

{a(x)=a0+a3x++an3xn31b(x)=a1+a4x+an2xn31c(x)=a2+a5x+an1xn31

我们定义 σ=n3 次单位根,那么有 ω=σ13,σ=ω3

于是只要我们先求出 a(x),b(x),c(x)σ0,,σn31 的值,就有

f(x)=a(x3)+xb(x3)+x2c(x3)f(ωi)=a(ω3i)+xb(ω3i)+x2c(ω3i)=a(σi)+xb(σi)+x2c(σi)

时间复杂度也是 O(nlogn)


同样地,我们可以把它推广到 NTT(也就是模 p 素数域)下面:

  • 卷积:n 放大到 3k(n3k),满足 3kp1,取 ω=gp13k
  • n 不能放大,np1,且 n 的质因子较小,可以根据 n 的最小质因子分治下去。

*单位根反演

难道他真的是天才?

[nv]=1nj=0n1(ωnv)j


证明是简单的:

由于

j=0n1xj=1xn1x(x1)

x=ωnv,则我们需要分两种情况讨论

  • x=1,则 1nj=0n1(ωnv)j=1,这种情况出现的条件是 nv
  • x1,则化简式的分母 1(ωnv)n=1(ωnn)v=11=0,而分子不为 0,所以原式 =0,这种情况出现的条件是 nv

得证。


单位根反演常常用来简化形如

v[nv]F(v)

的式子。

P.S.:这里的 ω 不一定是复数意义下的单位根,也可以是模 p 域的 n 次单位根。


有关 整除/同余 的题目 都可以 考虑 单位根反演!!!

有关 整除/同余 的题目 都可以 考虑 单位根反演!!!

有关 整除/同余 的题目 都可以 考虑 单位根反演!!!


上面扯了一堆看起来非常正确的东西,还是来练练题吧。


LOJ 6485. LJJ 学二项式定理 - 单位根反演

T 组数据,每次给定 n,s,a0,a1,a2,a3

i=0n(ni)siaimod4mod998244353

1T105,1n1018,1s,a0,a1,a2,a3108

煎蛋题。


看到 mod 想到 单位根反演

考虑答案为

t=03ati=0n(ni)si[4it]=14j=03(sω4j+1)nt=03atω4tj

结束!代码

Code
#include <bits/stdc++.h>
using namespace std;
#define int long long

const int H=998244353,g=3;
int n,s,a[4],w[4],ans=0,inv4,nw=0;

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int dec(int a,int b){return a<b?a-b+H:a-b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

void SOLVE(){
  cin>>n>>s>>a[0]>>a[1]>>a[2]>>a[3];
  ans=0;
  for(int j=0;j<4;j++){
	int nw=0;
	for(int t=0;t<4;t++) add(nw,mul(a[t],w[(4-t*j%4)%4]));
    add(ans,mul(qpow(adc(mul(s,w[j]),1),n),nw));
  }
  ans=mul(ans,inv4);
  cout<<ans<<'\n';
}

signed main(){
  /*2024.7.15 H_W_Y #6485. LJJ 学二项式定理 单位根反演*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  int _;cin>>_;
  w[1]=qpow(g,(H-1)/4),inv4=qpow(4);
  for(int i=2;i<=4;i++) w[i%4]=mul(w[i-1],w[1]);
  while(_--) SOLVE();
  return 0;
}

P5293 [HNOI2019] 白兔之舞 - 单位根反演

有一张顶点数为 (L+1)×n 的有向图。这张图的每个顶点由一个二元组(u,v)表示(0uL,1vn)
这张图不是简单图,对于任意两个顶点 (u1,v1)(u2,v2),如果 u1<u2,则从 (u1,v1)(u2,v2) 一共有 w[v1][v2] 条不同的边,如果 u1u2 则没有边。

白兔将在这张图上上演一支舞曲。白兔初始时位于该有向图的顶点 (0,x)

白兔将会跳若干步。每一步,白兔会从当前顶点沿任意一条出边跳到下一个顶点。白兔可以在任意时候停止跳舞(也可以没有跳就直接结束)。当到达第一维为 L 的顶点就不得不停止,因为该顶点没有出边。

假设白兔停止时,跳了 m 步,白兔会把这只舞曲给记录下来成为一个序列。序列的第 i 个元素为它第 i 步经过的边。

问题来了:给定正整数 ky1yn),对于每个 t0t<k),求有多少种舞曲(假设其长度为 m)满足 mmodk=t,且白兔最后停在了坐标第二维为 y 的顶点?

两支舞曲不同定义为它们的长度(m)不同或者存在某一步它们所走的边不同。

输出的结果对 p 取模。

p 为一个质数,108<p<2301n31xn1yn0w(i,j)<p1k65536kp1 的约数,1L108


考虑走 m 步的方案数,容易发现是

(Lm)Wmxy

于是,我们有

anst=m=0L(Lm)Wm[kmt]

直接单位根反演,再用二项式定理进行化简就可以得到

=1km=0L(Lm)Wmj=0k1(ωkmt)j=1kj=0k1ωktjm=0L(Lm)(ωkjW)m=1kj=0k1ωktj(I+ωkjW)L

我们令 Aj=(I+ωkjW)Laj=Aj[x][y],那么原式就是

1kj=0k1ωktjaj

发现这是一个 IDFT 形式,于是可以直接用 Bluestein 做了。


我们把 mk 展开成 (m+k2)(m2)(k2),代码较简单,直接 MTT 即可。代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ld long double
#define cp complex<ld>
#define ll long long
#define il inline

const int N=2.1e6+5,Bl=1<<15;
const ld PI=acos(-1);
int n,U[N],a[N],b[N],fac[N],ifac[N],H,prm[N],ct=0,g,ig,k,L,x,y,omg;
cp A[N],B[N],C[N],w[N];

il ll qpow(ll a,ll b=H-2){
  ll res=1;a%=H,b=(b%(H-1)+H-1)%(H-1);
  while(b){if(b&1) res=1ll*res*a%H;a=1ll*a*a%H,b>>=1;}
  return res;
}

il void FFT(cp *p,int len,int op){
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  cp nw,tmp,u;w[0]=cp(1,0);
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=cp(cos(2*PI/k),op*sin(2*PI/k));
    for(int i=1;i<=lst;i++) w[i]=w[i-1]*u;
	for(int i=0;i<len;i+=k)
	  for(int l=i;l<i+lst;l++)
	    tmp=w[l-i]*p[l+lst],p[l+lst]=p[l]-tmp,p[l]+=tmp;
  }
}

il void Init(int n){
  for(int i=0;i<n;i++) U[i]=(U[i>>1]>>1)|((i&1)?n>>1:0);
}

il void Mul(int *a,int *b,int n,int m){
  int P=1;
  for(P=1;P<n+m;P<<=1);
  Init(P);

  for(int i=0;i<n;i++) A[i]+=cp(a[i]/Bl,a[i]%Bl),C[i]+=cp(a[i]/Bl,-a[i]%Bl);
  for(int i=0;i<m;i++) B[i]+=cp(b[i]/Bl,b[i]%Bl);

  FFT(A,P,1),FFT(B,P,1),FFT(C,P,1);
  for(int i=0;i<P;i++) A[i]=A[i]*B[i],C[i]=C[i]*B[i];
  FFT(A,P,-1),FFT(C,P,-1);

  for(int i=0;i<P;i++){
	ll x=(ll)((A[i].real()+C[i].real())/P/2+0.49),y=(ll)(A[i].imag()/P+0.49),z=(ll)((C[i].real()-A[i].real())/P/2+0.49);
	a[i]=1ll*(1ll*x%H*Bl%H*Bl%H+1ll*y%H*Bl%H+z%H)%H;
  }
}

bool chk(int x){
  if(qpow(x,H-1)!=1) return false; 
  for(int i=1;i<=ct;i++)
    if(qpow(x,(H-1)/prm[i])==1) return false;
  return true;
}

int findrt(){
  for(int i=2;i<H;i++) if(chk(i)) return i;
  return -1;
}

void init(){
  int p=H-1;
  for(int i=2;i*i<=p;i++)
    if(p%i==0){
	  prm[++ct]=i;
	  while(p%i==0) p/=i;
	}
  if(p!=1) prm[++ct]=p;

  g=findrt(),ig=qpow(g);
}

struct Mat{
  int a[4][4];
  void init(){memset(a,0,sizeof(a));}
}S;

Mat operator *(Mat a,Mat b){
  Mat c;c.init();
  for(int k=1;k<=n;k++)
    for(int i=1;i<=n;i++)
      for(int j=1;j<=n;j++)
	    (c.a[i][j]+=1ll*a.a[i][k]*b.a[k][j]%H)%=H;
  return c;
}

Mat operator *(Mat a,int b){
  for(int i=1;i<4;i++) for(int j=0;j<4;j++) a.a[i][j]=1ll*a.a[i][j]*b%H;
  return a;
}

Mat ins(Mat a){
  for(int i=1;i<=n;i++)
    (a.a[i][i]+=1)%=H;
  return a;
}

Mat Qpow(Mat a,int b){
  Mat c;c.init();
  for(int i=1;i<=n;i++) c.a[i][i]=1;
  while(b){if(b&1) c=c*a;a=a*a,b>>=1;}
  return c;
}

int bn(int x){return 1ll*x*(x-1)/2%(H-1);}

void SOLVE(){
  omg=qpow(g,(H-1)/k);
  int nw=1;
  for(int i=0;i<k;i++,nw=1ll*nw*omg%H)
	a[i]=1ll*Qpow(ins(S*nw),L).a[x][y]*qpow(omg,bn(i))%H;

  k<<=1;
  for(int i=0;i<k;i++) b[i]=qpow(omg,H-1-bn(i));
  reverse(b,b+k);
  Mul(a,b,k,k);
  reverse(a,a+k);

  int I=qpow(k>>1);
  for(int i=0;i<(k>>1);i++) cout<<1ll*I*a[i]%H*qpow(omg,bn(i))%H<<'\n';
}

signed main(){
  /*2024.7.15 H_W_Y P5293 [HNOI2019] 白兔之舞 单位根反演 + Bluestein*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>k>>L>>x>>y>>H;
  for(int i=1;i<=n;i++)
    for(int j=1;j<=n;j++)
	  cin>>S.a[i][j];
  init();
  SOLVE();
  return 0;
}

这份代码 是展开成上面 Bluestein 的第一个形式,但是挂了,求调/kk


P5591 小猪佩奇学数学 - 单位根反演

给定 n,p,k,询问

i=0n(ni)×pi×ikmod998244353

1n,p<998244353,k{2w|0w20}

居然切掉了/jk


很明显不能整除分块,于是我们考虑把下取整拆开

i=0n(ni)×pi×ik=i=0n(ni)×pi×iimodkk=1k[i=0n(ni)ipii=0n(ni)(imodk)pi]

对于前者,我们考虑吸收公式

i=0n(ni)ipi=i=0n(n1i1)npi=npi=0n1(n1i)pi=np(p+1)n1

而对于后者,我们直接单位根反演

i=0n(ni)(imodk)pi=t=0k1ti=0n(ni)pi[kit]=1kt=0k1ti=0n(ni)pij=0k1(ωkit)j=1kt=0k1tj=0k1ωktji=0n(ni)(pωkj)i=1kt=0k1tj=0k1ωktj(pωkj+1)n

仔细一看,后面的 j 中不就是 (pωkj+1)n 的 IDFT 吗?

由于保证了 k=2ω,所以可以直接 NTT。

于是就做完了,时间复杂度 O(klogk)。还跑得蛮快的。代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ll long long 
#define ull unsigned long long

const int N=2.1e6+5,H=998244353;
int n,p,m,invm,omg,ans=0,U[N];
ull F[N],w[N];

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int dec(int a,int b){return a<b?a-b+H:a-b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}
void del(int &a,int b){a=dec(a,b);}

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
  return res;
}

const int g=3,ig=qpow(g);

void NTT(ull *p,int len,int op){
  for(int i=0;i<len;i++) U[i]=(U[i>>1]>>1)|((i&1)?len>>1:0);
  for(int i=0;i<len;i++) if(i<U[i]) swap(p[i],p[U[i]]);
  ull u,tmp;w[0]=1;
  for(int k=2,lst=1;k<=len;lst=k,k<<=1){
	u=qpow(op==1?g:ig,(H-1)/k);
    for(int i=1;i<=lst;i++) w[i]=w[i-1]*u%H;
	for(int i=0;i<len;i+=k)
	  for(int l=i;l<i+lst;l++)
	    tmp=w[l-i]*p[l+lst]%H,p[l+lst]=p[l]-tmp+H,p[l]+=tmp;
	if(k==(1<<16)) for(int i=0;i<len;i++) p[i]%=H;
  }
  for(int i=0;i<len;i++) p[i]%=H;
}

int main(){
  /*204.7.15 H_W_Y P5591 小猪佩奇学数学 单位根反演*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  cin>>n>>p>>m;invm=qpow(m);
  omg=qpow(g,(H-1)/m);
  ans=mul(mul(n,p),qpow(adc(p,1),n-1));
  for(int i=0,nw=1;i<m;i++,nw=mul(nw,omg)) F[i]=qpow(adc(mul(p,nw),1),n);
  NTT(F,m,-1);
  for(int i=0;i<m;i++) del(ans,mul(mul(i,invm),F[i]));
  ans=mul(ans,invm);
  cout<<ans<<'\n';
  return 0;
}

CQYC 7.13 联考 - B.人じゃありませんから(jinrui)- 单位根反演

奈亚子很喜欢集合。

有一天真寻给了她一个集合 S={1,2,,n},她想知道有多少子集 TS,满足 gcd(tTt,n)=1

通俗地说,求集合 {1,2,,n} 中,有多少子集的元素之和与 n 互素。

由于奈亚子设立了方便剧情发展的结界,所以你可以对 109+7 取模。

奈亚子因为可以生命时间加速,一瞬间就做出来了,而你只有 2s 的时间,你还需要处理 T 组数据。

T10,n2×1015,保证 n 不为 109+7 的倍数。

乎✌的胜利!知乎


m=n(n+1)2,那么问题即求

i=1m[gcd(i,n)=1]f(i)

其中 f(i) 表示和为 i 的子集个数。

而我们发现 f(i) 的生成函数其实是 F(x)=i=1n(1+xi),也就有 f(i)=[xi]F(x)

我们对这个式子进行莫比乌斯反演

i=1m[gcd(i,n)=1]f(i)=i=1mdn,diμ(d)f(i)=dnμ(d)i=1m[di]f(i)

看到 [di],我们想打单位根反演,于是有

i=1m[di]f(i)=1di=1mj=0d1(ωdi)jf(i)=1dj=0d1i=1m(ωdj)i[xi]F(x)=1dj=0d1F(ωdj)

现在式子就变成了求

dnμ(d)di=0d1F(ωdi)

那么我们现在想去化简后面的部分,也就是

i=0d1F(ωdi)=i=0d1(j=1n(1+ωdij))

发现上面的 ij 是一个经典问题,也就是设 gcd(i,d)=k,则 ωdij 会循环 k 次,每次循环节长度是 dk,那么可以得到

i=0d1(j=1n(1+ωdij))=kd(i=0d1[gcd(i,d)=k])(i=1dk(1+ωdik))nkd=kdφ(dk)(i=1dk(1+ωdik))nkd

然后考虑结论

i=1n(xωni)=xn1

证明就是考虑 xn1n 个根是 ωni,于是我们又可以将原式化简

kdφ(dk)(i=1dk(1+ωdik))nkd=kdφ(dk)((1)dki=1dk(1ωdik))nkd=kdφ(dk)(1(1)dk)nkd=kdφ(k)(1(1)k)nk=kd,2kφ(k)2nk

于是答案就是

dnμ(d)dkd,2kφ(k)2nk=kn,2kφ(k)2nkdn,kdμ(d)d

我们考虑继续化简后面的那一块。

注意到一些数论函数的式子

μ(ij)=μ(i)μ(j)[gcd(i,j)=1]φ(ij)=φ(i)φ(j)gcd(i,j)φ(gcd(i,j))φ(i)=φ(ij)φ(gcd(i,j))φ(j)gcd(i,j)dnμ(d)d=φ(n)n

都是可以用定义式或者计算式简单证明出来的,于是我们可以继续推导

dn,kdμ(d)d=dnkμ(kd)kd=μ(k)kdnkμ(d)[gcd(k,d)=1]d=μ(k)kdnkgcd(k,nk)μ(d)d=μ(k)kkgcd(k,nk)φ(nkgcd(k,nk))n

x=kgcd(k,gcd(k,nk)),于是

φ(nx)=φ(n)φ(gcd(nx,x))φ(x)gcd(nx,x)

由于 gcd(kgcd(k,nk),nkgcd(k,nk))=1,所以

φ(nx)=φ(n)φ(x)

而我们还可以继续化简 φ(x),也就是

φ(x)=φ(kgcd(k,nk))=φ(k)φ(gcd(k,nk))gcd(k,gcd(k,nk))φ(gcd(k,gcd(k,nk)))

稍微化简一下就是

φ(x)=φ(k)gcd(k,nk)

于是式子就可以继续化简

μ(k)kkgcd(k,nk)φ(nkgcd(k,nk))n=μ(k)kkgcd(k,nk)φ(n)nφ(k)gcd(k,nk)=φ(n)nμ(k)φ(k)

那么

ans=kn,2kφ(k)2nkdn,kdμ(d)d=φ(n)nkn,2kφ(k)2nkμ(k)φ(k)=φ(n)nkn,2kμ(k)2nk

就做完了!实现非常简单。

注意特判 n=1 时答案是 1代码

Code
#include <bits/stdc++.h>
using namespace std;
#define int long long

const int N=20,H=1e9+7;
int cnt=0,ans=0,n,a[N],p,phi;

int qpow(int a,int b=H-2){
  int res=1;
  while(b){if(b&1) res=res*a%H;a=a*a%H,b>>=1;}
  return res;
}

void dfs(int id,int s,int v){
  if(id>cnt) return ans=(ans+v*qpow(2,(n/s))%H)%H,void();
  dfs(id+1,s,v),dfs(id+1,s*a[id],H-v);
}

void SOLVE(){
  cin>>n;cnt=ans=0;
  p=phi=n;
  if(n<=2) return cout<<n<<'\n',void();
  if(!(p&1)) phi/=2;
  while(!(p&1)) p/=2;
  for(int i=3;i*i<=p;i++)
    if(!(p%i)){
	  phi=phi/i*(i-1),a[++cnt]=i;
	  while(!(p%i)) p/=i;
	}
  if(p!=1) phi=phi/p*(p-1),a[++cnt]=p;
  dfs(1,1,1);
  cout<<(phi%H*qpow(n%H)%H*ans%H)<<'\n';
}

signed main(){
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  int _;cin>>_;
  while(_--) SOLVE();
  return 0;
}

P10084 [GDKOI2024 提高组] 计算 - gcd + 单位根反演

定义 F(x,a,b)=gcd(xa1,xb1)+1,x>0

特别的,如果 a=0b=0F(x,a,b)=0

现在给出五个非负整数 m,a,b,c,d

L=F(m,a,b)+1R=F(m,c,d)

问集合 {L,L+1,L+2,,R2,R1,R} 有多少个子集和是 m 的倍数。

由于答案可能很大,你只需要输出方案数对 998244353 取模后的结果就可以了。

0<m107,L,R1018,a,b,c,d103,T104

如果做了上一道题,那么这道题就跟坎瓜切菜一样。


考虑 L,R 究竟是什么?

对于 gcd,我们不仅要考虑 辗转相除法,还要去思考 辗转相减法

在这道题中,用 辗转相减 可以得到

gcd(xa1,xb1)=gcd(xb1,xb(xab1))

然而,显然 xb1xb,所以上式就变成了

gcd(xb1,xab1)

发现系数也是更相减损!这也是一个 gcd 的过程,到最后

gcd(x01,xp1)=xp1

于是答案其实是

xgcd(a,b)

打表应该也能看出来


那么现在我们知道了 L=mgcd(a,b)+1,R=mgcd(c,d),容易得知 mRL+1,这一点会在之后的推到中有用。

于是类似于上一个题,我们就有了一下的推导。

首先列出答案式子

ans=i=0sum[mi]f(i)

其中 sum=(L+R)(RL+1)2,f(i) 表示多少个集合的和是 i

容易知道 f(i) 的生成函数是

F(x)=i=LR(1+xi)

我们对原式做 单位根反演,可以得到

ans=i=0sum[mi]f(i)=1mi=0sumj=0m1(ωmi)jf(i)=1mi=0m1F(ωmi)=1mi=0m1(j=LR(1+ωmij))=1mi=0m1(j=0m1(1+ωmij))RL+1m=1mdm(i=0m1[gcd(i,m)=d])(j=0md(1+ωmdj))dRL+1m=1mdmφ(d)(1(1)d)RL+1d=1mdm,2dφ(d)2RL+1d

以上推导过程较为简略,但基本上和上一道题一样,看不懂的地方可以去看上一道题。


于是这样我们就得到了答案,直接枚举因数即可。

注意空集也要算,所以 L>R 的情况答案是 1代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ll long long

const int N=105,H=998244353;
int prm[N],c[N],ct=0,m,A,B,C,D,ans=0,x;
ll L,R;

int adc(int a,int b){return a+b>=H?a+b-H:a+b;}
int mul(int a,int b){return 1ll*a*b%H;}
void add(int &a,int b){a=adc(a,b);}

int gcd(int a,int b){return !b?a:gcd(b,a%b);}

ll qpow(ll a,ll b=H-2,ll op=0){
  if(!op){
	ll res=1;a%=H,b%=(H-1);
	while(b){if(b&1) res=mul(res,a);a=mul(a,a),b>>=1;}
	return res;
  }
  ll res=1;
  while(b){if(b&1) res=res*a;a=a*a,b>>=1;}
  return res;
}

ll F(int x,int a,int b){
  if(!a||!b) return 0;
  return qpow(x,gcd(a,b),1);
}

void dfs(int id,int phi,int v){
  if(id==ct+1) return add(ans,mul(phi,qpow(2,(R-L+1)/v)));
  dfs(id+1,phi,v);
  for(int i=1,nw=1;i<=c[id];i++,nw=mul(nw,prm[id]))
	dfs(id+1,mul(phi,mul(nw,prm[id]-1)),mul(v,mul(nw,prm[id])));
}

void SOLVE(){
  cin>>m>>A>>B>>C>>D;x=m,ct=ans=0;
  L=F(m,A,B)+1,R=F(m,C,D);
  if(L>R) return cout<<"1\n",void();

  while(x%2==0) x/=2;
  for(int i=3;i*i<=x;i++){
	if(x%i==0){
	  prm[++ct]=i,c[ct]=0;
	  while(x%i==0) ++c[ct],x/=i;
	}
  }
  if(x!=1) prm[++ct]=x,c[ct]=1;
  dfs(1,1,1);

  cout<<mul(qpow(m),ans)<<'\n';
}

int main(){
  /*2024.7.19 H_W_Y P10084 [GDKOI2024 提高组] 计算 单位根反演*/
  ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
  int _;cin>>_;
  while(_--) SOLVE();
  return 0;
}

P3746 [六省联考 2017] 组合数问题 - 循环卷积

给定 n,p,r,k(1n109,0r<k50,2p2301),求

imodk=r(nki)modp

不知道为什么要用单位根反演推。


直接变成 (1+x)nk 的循环卷积即可。代码

Code
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define poly vector<int>

int p,k,n,r;

poly operator *(poly a,poly b){
  poly c(k);
  for(int i=0;i<k;i++)
    for(int j=0;j<k;j++)
	  (c[(i+j)%k]+=1ll*a[i]*b[j]%p)%=p;
  return c;
}

int main(){
  /*2024.7.15 H_W_Y P3746 [六省联考 2017] 组合数问题 循环卷积*/
  cin>>n>>p>>k>>r;
  poly a(k),ans(k);
  if(k==1) a[0]=2%p;
  else a[0]=a[1]=1;
  ans[0]=1;

  ll b=1ll*n*k;
  while(b){if(b&1) ans=ans*a;a=a*a,b>>=1;}
  cout<<ans[r]<<'\n';
  return 0;
}

本文作者:H_W_Y

本文链接:https://www.cnblogs.com/H-W-Y/p/18317563/poly1

版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。

posted @   H_W_Y  阅读(35)  评论(1编辑  收藏  举报
点击右上角即可分享
微信分享提示
评论
收藏
关注
推荐
深色
回顶
收起