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(1-x) = - \sum_{i=1}^{\infty} \frac {x^i} i \]

证明两边同时求导即可。


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


*MTT

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


其核心是把多项式系数 \(F_i\) 拆成 \(A_ik+B_i\) 的形式，从而解决 值域过大 和 精度过高 的问题。

例如在 P4245 【模板】任意模数多项式乘法 中，乘出来的值高达 \(10^9 \times 10^9 \times 10^5 = 10^{23}\)，我们尝试取 \(k=2^{15}\)，那么就可以把 \(A,B\) 的最终值域优化到 \(10^{14}\)，非常优秀。

那么这道题就变成了求

\[\begin{aligned} & (A_0 k + A_1)(B_0 k + B_1) \\ = & A_0 B_0 k^2 + (A_0B_1+A_1B_0)k+A_1B_1 \end{aligned} \]

直接做可以做到用 \(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 中的三步变两步，将其转化成复数多项式。

也就是我们可以设

\[\begin{aligned} P & = A_0 + A_1 i \\ Q & = B_0 + B_1 i \\ P' & = A_0 - A_1 i \end{aligned} \]

那么就可以得到

\[\begin{aligned} P*Q & = A_0B_0 -A_1B_1 + (A_0B_1+A_1B_0) \cdot i \\ P'*Q & = A_0B_0 + A_1B_1 +(A_0B_1-A_1B_0) \cdot i \end{aligned} \]

于是我们需要的东西都可以从上面两个东西做加减出来了，容易发现这样只需要 \(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\) 个数 \(q_1,q_2, \dots q_n\)，定义

\[F_j~=~\sum_{i = 1}^{j - 1} \frac{q_i \times q_j}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i \times q_j}{(i - j)^2},E_i~=~\frac{F_i}{q_i} \]

对 \(1 \leq i \leq n\)，求 \(E_i\) 的值。

\(1 \leq n \leq 10^5\)，\(0 < q_i < 10^9\)。

典题先来一个。


首先推一下式子

\[E_j = \sum_{i \lt j} \frac {q_i} {(i-j)^2} - \sum_{i \gt j} \frac {q_i} {(i-j)^2} \]

我们设

\[d_k = \begin{cases} \frac 1 {k^2} & k \gt 0 \\ 0 & k=0 \\ - \frac 1 {k^2} & k \lt 0 \end{cases} \]

于是我们有

\[E_j = \sum_{i=0}^n q_i d_{j-i} \]

但是这不是标准卷积形式，所以我们考虑 对称（镜像）。

令 \(E'_{j+n}=E_j,D'_{k+n}=D_k\)，于是就有了以下式子

\[E'_{j+n} = \sum_{i+k=j+n} q_i D'_k \]

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

注意精度问题，我们需要给 \(D\) 乘上 \(10^{11}\)。代码。

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 \times C\) 的棋盘，有一些”超级车“，问：有多少个格子未被攻击。

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

\(R,C,N \le 10^5\)。

小容斥一手。


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

于是我们设

\[\begin{cases} R:被同行的车攻击的格子集合 \\ C:被同列的车攻击的格子集合 \\ D:被同对角线的车攻击的格子集合 \end{cases} \]

那么其实答案是

\[|R \cup C \cup D |= |R| + |C| + |D| - |R \cap C| - |R \cap D| - |C \cap D| + |R \cap C \cap D| \]

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

于是其实只用讨论如何求 \(|R \cap C \cap D|\) 了。


不妨设

\[\begin{cases} r_1,r_2,\cdots,r_n 行有车 \\ c_1,c_2,\cdots,c_m 列有车 \\ d_1,d_2,\cdots,d_l 对角线有车 \end{cases} \]

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

\[r_i + d_k = c_j \]

这种和式不就是 多项式卷积 形式？！

于是我们令

\[R(x) = \sum x^{r_i},D(x) = \sum x^{d_k} \]

那么答案就是

\[|R \cap C \cap D| = \sum_{j} [x^{c_j}] R(x)D(x) \]

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

时间复杂度 \(\mathcal O(n \log n)\)。代码。

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个整数\(A_1,A_2,\dots,A_N\),对所有的\(S\)，求满足$A_i+A_j+A_k = S \(且\)i < j < k\(的\)(i,j,k)$对数.

\(N \le 40000, A_i \le 20000\)


首先来考虑如果没有 \(i,j,k\) 的限制怎么办？

容易发现其实就是

\[[x^S]\left(\sum_{i=1}^N x^{A_i} \right )^3 \]

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

这个东西其实是

\[\left ( \sum_{i=1}^N x^{A_i} \right)^3= \sum_{i=1}^N x^{3A_i} +3 \sum_{i \neq j} x^{2A_i+A_j} + 6 \sum_{i \lt j \lt k} x^{A_i+A_j+A_k} \]

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

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

\[\left ( \sum_{i=1}^N x^{2A_i} \right) \left( \sum_{i=1}^N x^{A_i} \right) \]

再来分讨一下它是什么？

容易发现其实就是

\[\left ( \sum_{i=1}^N x^{2A_i} \right) \left( \sum_{i=1}^N x^{A_i} \right) = \sum_{i=1}^N x^{3A_i} + \sum_{i \neq j} x^{2A_i+A_j} \]

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

于是就做完了，时间复杂度 \(\mathcal O(V \log V)\)。代码。

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 \times N\) 的网格图。一些点会落在网格图上，记 \(C_{xy}\) 为坐标是 \((x, y)\) 的点的数量 (\(0 \le C_{xy} \le 9\))。现要求对于所有 成对的点之间的欧几里得距离，统计：

  • 全部距离依次平方后，排序并去重后的结果 \(d_1, d_2, \ldots, d_m\)；
  • 上述序列去重前 \(d_i\) 重复出现的次数 \(c_i\)；
  • 所有距离的平均距离 \(D_{ave}\)。

\(1 \le N \le 1024\)，\(0 \le C_{xy} \le 9\)。

应该是紫吧？


先写一个式子出来

\[D_{x,y} = \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} C_{i,j} C_{i+x,j+y} \]

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

\[\phi (x,y) = x \times 2N+y \]

同时令 \(M=2N^2\)，那么就可以得到

\[D_{\phi(x,y)}=\sum_{i=0}^{N-1} \sum_{j=0}^{N-1} C_{\phi(i,j)} C_{\phi(i,j)+\phi(x,y)} \]

于是有

\[D_z = \sum_{u=0}^M C_uC_{u+z} \]

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

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

\[\begin{cases} C'_k = C_{M-k}\\ D'_k = D_{M-k}\\ \end{cases} \]

于是就得到了

\[D'_k = \sum_{i=0}^k C_i C'_{k-i} \]

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

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

时间复杂度 \(\mathcal O(M \log M)\)。代码。

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\)（分别用 \(pS\)、\(pH\)、\(pC\) 和 \(pD\) 表示）。没有其它类型的牌。

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

390770

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

\(b \le 50000,c \le 10000\)。

唐唐板子题。


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

注意每组数据之后要输出一个 \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\) 组的方案数。

\(n\le 10^9\)，\(k<2^{15}\)。


首先，先写一个 dp，即设 \(f_{i,j}\) 表示前 \(i\) 个数，分成 \(j\) 组方案数。

于是我们有

\[f_{i,j} = f_{i-1,j} + f_{i-1,j-1} + f_{i-2,j-1} \]

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

\[F_i(x) = (x+1) F_{i-1}(x)+xF_{i-2}(x) \]

我的天，矩阵，快速，幂？

于是构造矩阵

\[A=\left [\begin{matrix} x+1 & 1 \\ x & 0 \end{matrix} \right] \]

那么答案其实就是 \(A^n[0][0]\)。

做完了，直接暴力 NTT 即可，常数大，时间复杂度 \(\mathcal O(k \log n \log k)\)。

笔者本以为它过不了——但是它过了！需要开 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。

考虑能不能从 \(f_a\) 和 \(f_b\) 转移到 \(f_{a+b}\)？

发现是可以的

\[f_{a+b,i} = \sum_{j=0}^i f_{a,j} f_{b,j-i} + \sum_{j=0}^{i-1} f_{a-1,j} f_{b-1,i-j-1} \]

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

\[F_{a+b}(x) = F_a(x)F_b(x)+xF_{a-1}(x)F_{b-1}(x) \]

于是我们尝试从 \(n \to 2n\)，也就是

\[\begin{cases} F_{2n} (x)= F_n(x)^2+xF_{n-1}(x)^2 \\ F_{2n-1}(x) = F_n(x)F_{n-1}(x)+x F_{n-1}(x) F_{n-2}(x) \\ F_{2n-2}(x) = F_{n-1}(x)^2 + F_{n-2}(x)^2 \end{cases} \]

那么每次我们只需要维护 \(F_n(x),F_{n-1}(x),F_{n-2}(x)\) 就可以通过上述转移转移到 \(F_{2n}(x),F_{2n-1}(x),F_{2n-2}(x)\)。

这就是 倍增 NTT 了，直接做还是 \(\mathcal O(k \log n \log k)\)。没有实现（应该和上述代码类似），咕。


另外，这里还有一种 \(\mathcal O(k \log k)\) 的高级做法。

观察上面的式子

\[F_i(x) = (x+1) F_{i-1}(x)+xF_{i-2}(x) \]

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

也就是设

\[G(z) = F_0 + F_1z+ F_2 z^2 + \cdots \]

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

\[G(z) = (x+1)G(z)z+xG(z)z^2+1 \]

转化出来可以得到

\[G(z) = \frac {1} {1-(1+x)z-xz^2} \]

我们需要求出 \(F_n\)，于是想把这个式子展开，展开成

\[\frac u {1-Az} + \frac {v} {1-Bz} \]

的形式。


于是直接列方程

\[\begin{cases} A + B = 1+x \\ AB=-x \end{cases} 解之得： \begin{cases} A=& \frac {(1+x)+\sqrt{(1+x)^2+4x}} 2 \\ B = & \frac {(1+x)-\sqrt{(1+x)^2+4x}} 2 \end{cases} \]

带入 \(u,v\) 的方程

\[\begin{cases} u+v=1 \\ uB+vA=0 \end{cases} 解之得： \begin{cases} u = & \frac A{A-B} \\ v= & \frac {-B}{A-B} \end{cases} \]

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

\[\begin{aligned} G(z) = & u \cdot \sum_{n=0}^\infty A^n z^n + v \cdot \sum_{n=0}^{\infty} B^n z^n \\ [z^n] G(z) = & u A^n+ vB^n = \frac {A^{n+1}-B^{n+1}} {A-B} \end{aligned} \]

注意到 \(B^{n+1} \equiv 0 \pmod {x^{n+1}}\)，所以其实答案式子是

\[\frac {A^{n+1}} {A-B} \]

最后的式子就可以以 \(\mathcal O(k \log k)\) 的时间求得了。代码。

这是一个 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 在字符串匹配中的应用

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

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

\(1 \le m \le n \le 3 \times 10^5\)。

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


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

那怎么办呢？

考虑定义配对函数

\[C_{x,y}= [A_x-B_y]^2A_xB_y \]

当前字母是通配符时 \(A_x =0\)，容易发现，当且仅当 \(x,y\) 匹配时 \(C\) 为 \(0\)。

考虑设 \(P_k\) 表示 \(k\) 位置为结尾的匹配情况，于是 \(k\) 位置时合法的，当且仅当

\[P_k = \sum_{i=0}^{m-1} [A_i-B_{k-m+i+1}]^2A_iB_{k-m+i+1} =0 \]

这并不好看。考虑按照之前的套路，我们令 \(S_i = A_{m-1-i}\)，也就可以得到

\[P_k = \sum_{i=0}^{m-1} [S_i-B_{k-i}]^2 S_iB_{k-i} \]

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

\[P_k = \sum_{i=0}^k {S_i}^3B_{k-i} + \sum_{i=0}^k S_i{B_{k-i}}^3 - 2\sum_{i=0}^k {S_i}^2 {B_{k-1}}^2 \]

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

于是就做完了，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\) 和两个只包含 \(\texttt{AGCT}\) 四种字符的基因串 \(S\) 和 \(T\)。现在你要找出在下列规则中 \(T\) 在 \(S\) 中出现了几次。

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

即对于所有的 \(j \in[1,|T|]\)，都存在一个 \(p \in [1,|S|]\) 使得 \(|(i+j-1)-p| \leq k\) 且 \(S_p=T_j\) 。

例如 \(k=1\) 时，\(\texttt{ACAT}\) 出现在 \(\texttt{AGCAATTCAT}\) 的 \(2\) 号， \(3\) 号和 \(6\) 号位置。 (编号从 \(1\) 开始。)

\(1≤n≤m≤2*10^5\ ,\ 0≤k≤2*10^5\)

居然切了。


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

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

  • \(A\)：若 \(A_i = 1\) 则表示 \(i\) 位置合法，即存在一个和 \(i\) 距离不超过 \(k\) 的碱基；反之则 \(A_i = 0\)。
  • \(B\)：若 \(B_i = 1\) 则表示 \(T_{m-1-i}\) 是当前的碱基，反之则不是。

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

\[P_k = \sum_{i=0}^k B_i A_{k-i} =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 - 拉插

求 \((\sum_{i=1}^ni^k) \bmod (10^9+7)\)。

\(1 \le n \le 10 ^ 9, 0 \le k \le 10 ^ 6\)。

有些人拉插喜欢把 \(x_i\) 当成 \(y_i\) 算。


首先，为什么这是一个 \(k+1\) 次多项式？

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

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

由于我们把这个式子差分之后是 \(1^k,2^k,\cdots,n^k\) 是 \(k\) 次多项式，所以原式是一个 \(k+1\) 次多项式。


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

\[\sum_{i=1}^{k+2} y_i \prod_{j \neq i} \frac {n-j} {i-j} \]

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

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

对于一个整数序列 \(\{a_1, a_2, \ldots, a_n\}\)，我们定义它的变换为 \(\{b_1, b_2, \ldots, b_n\}\)，其中 \(b_i = a_1 | a_2 | \ldots | a_i\)，其中 \(|\) 表示二进制按位或运算。

现在求有多少个长为 \(n\) 的序列 \(a\)，满足 \(1\leq a_i \leq 2^k - 1\)，使得它的变换 \(b\) 是严格单调递增的，对 \(10^9+7\) 取模。

\(1\leq n \leq 10^{18}\)，\(1\leq k \leq 3 \times 10^4\)。

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


首先，容易发现 \(n \gt k\) 的答案是 \(0\)，于是接下来我们默认 \(n \le k\) 了。

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

那么可以得到

\[f_{i,j} = \sum_{p=1}^j f_{i-1,j-p} 2^{j-p} \binom j p \]

看起来就非常卷积。

于是考虑倍增 NTT，考虑知道 \(f_a,f_b\) 时如何推出 \(f_{a+b}\)，容易列出式子

\[\begin{aligned} f_{a+b,j} & = \sum_{p=1}^j f_{a,p} f_{b,j-p} \binom j p 2^{bp} \\ \frac {f_{a+b,j}} {j!} & = \sum_{p=0}^j \frac {f_{a,p} 2^{bp}} {p!} \frac {f_{b,j-p}} {(j-p)!} \end{aligned} \]

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

时间复杂度 \(\mathcal O(k \log k\log n)\)。代码。

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\)。并且由于这个手环是一个圆，可以以任意的角度旋转它，但是由于上面装饰物的方向是固定的，所以手环不能翻转。需要在经过亮度改造和旋转之后，使得两个手环的差异值最小。

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

\[\sum_{i=1}^{n} (x_i-y_i)^2 \]

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

\(1 \le n \le 50000\), \(1 \le a_i \le m \le 100\)。

居然自己做出来了 /jk


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

\[P_i = \sum_{j=0}^{n-1} (a_j - b_{j+i})^2 \]

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

于是我们设

\[P'_i = P_{2n-1-i},b'_i = b_{2n-1-i} \]

那么

\[{P'}_i = \sum_{j=0}^{n-1} {a_j}^2 + \sum_{j=0}^{n-1}{b'}_{i-j}^2 - 2 \sum_{j=0}^{i} a_j {b'}_{i-j} \]

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

而假设变化量是 \(c\)，我们直接在 \(b'\) 的每一项加上一个 \(c\)，拆下括号发现是可以 \(\mathcal 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 \times 2 ^{21} + 1\) ) 取模即可。

\(n \le 130000\)。

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


考虑经典结论，若 \(A\) 中包含若干个 \(B\)，则 \(\hat A(x) = \exp \left(\hat B (x) \right)\)。

其中 \(\hat A(x)\) 表示 \(A\) 的指数生成函数，\(\hat B(x)\) 同理。

于是设 \(\hat A(x)=\sum_{n=0}^{\infty} \frac {2^{\binom n 2}} {n!} x^n\)，即带标号的图个数，于是答案 \(\hat B(x) = \ln \left(\hat A(x) \right)\)，结束！

时间复杂度 \(\mathcal O(n \log n)\)。代码。

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

这个背包最多可以装 \(10^5\) 大小的东西

付公主有 \(n\) 种商品，她要准备出摊了

每种商品体积为 \(v_i\)，都有无限件

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

\(1\le n,m \le 10^5\)，\(1\le v_i \le m\)。

\(1.17s\) 略胜一筹。


考虑设 \(a_i\) 表示体积为 \(i\) 的物品种类数，于是我们有生成函数

\[\prod_{1 \le i \le m} (1+x^i+x^{2i}+\cdots)^a_i = \prod_{1 \le i \le m} \left( \frac 1 {1-x^i} \right)^{a_i} \]

不知道怎么做？经典操作：先去 Ln，再去 Exp，可以得到

\[\begin{aligned} & \prod_{1 \le i \le m} \left( \frac 1 {1-x^i} \right)^{a_i} \\ = & \exp \left( - \sum_{i=1}^m a_i \ln \left( 1-x^i \right) \right) \\ = & \exp \left( \sum_{i=1}^m a_i \sum_{j=1}^{\left \lfloor \frac m i \right \rfloor} \frac {x^{ij}} j \right) \end{aligned} \]

后面的式子是调和级数形式，只有 \(O(n \log n)\) 项，可以暴力求和。

于是最后再做一次 \(\exp\) 就可以了，时间复杂度 \(\mathcal O(m \log m)\)。代码。

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\times 17 \times 2^{23}+1)\) 一个质数意义下，给出了 \(m\) 条指令，每一条指令都是下列两种之一：

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

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

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

\(1 \le n \le 10^5,1 \le m \le 60000\)。

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


我们设初始的序列为 \(\left < x_1,x_2,\cdots,x_n \right>\)，那么第 \(i\) 次操作后，\(x_j\) 一定形如

\[\frac {a_i x_j + b_i} {c_i x_j +d_i} \]

而这里的 \(a_i,b_i,c_i,d_i\) 是可以线性求出来的。

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

\[\frac 1 {x_j + t_i} \times k_i + e_i \]

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

\[F(t) = \sum_{i=1}^n \frac 1 {x_i+t} \]

的 \(t=t_i\) 时的值，那么这就是一个多点求值了。


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

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

\[\prod_{i=1}^n (x_i +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 有一棵树。这棵树上只有点 \(a\) 和 \(b\) 是黑色，其它的点都是白色。

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

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

\(a,b\le n\le 234,567\)。

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


考虑不在 \(a \to b\) 路径上的子树，那么每一次一定是最先染根，这就是一个满足树形拓扑序的排列问题。

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

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


于是按照这个思路想下去，我们关心地其实就是 \(a \to b\) 路径上点染色的顺序。

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

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

先给环上点标号，\(s\) 为 \(0\)，其余依次从 \(1 \sim k\)，其中 \(k\) 表示环上的点的个数。


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

怎么搞都是 \(n^2\) 的 dp，怎么办呢？

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

发现是可以的，我们枚举基环树上的一条边 \(r \to r+1\)，断掉它表示最后一个点要么是 \(r\) 要么是 \(r+1\)，于是最后需要乘上 \(\frac 1 2\)。

我们再钦定每个点的权值是 \(a_i\) 表示环上第 \(i\) 个点，其非环上的子树大小之和加一；\(S_i\) 为 \(a_i\) 的前缀和。

那么答案就是

\[\frac 1 2 C \left( \prod_{i=0}^k \frac 1 {\prod_{j=1}^i (S_i - S_{j-1}) \prod_{j=i+1}^k (S_j - S_i)} \right) \]

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

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

\[\prod_{j=1}^i (S_i - S_{j-1}) \prod_{j=i+1}^k(S_j-S_i) \]

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

\[F_i(x) = \frac {\prod_{j=0}^k (x-S_j)} {x-S_i} \]

在 \(x \to S_i\) 时的值再乘上一个 \((-1)^{k-i}\)。

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

\[F(x) = \left( \prod_{j=0}^k (x-S_j) \right)' \]

在 \(x=S_i\) 处的值。


由于我们要同时求每个 \(x=S_i\) 的值，所以这是需要用多项式多点求值完成的。

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

\[ans = \frac 1 2 C \left( \sum_{i=0}^k \frac 1 {F(S_i) (-1)^{k-i}} \right) \]

时间复杂度 \(\mathcal O(n \log^2 n)\)。代码。

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) = \sum_{i=0}^n a_i x^{i \bmod k} \]

从 \(n\) 次多项式变成 \(k\) 次多项式。

然而这两个多项式满足

\[A(\omega_k^0)\cdots A(\omega_k^{k-1}) = B(\omega_{k}^0) \cdots B(\omega_{k}^{k-1}) \]

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


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

\[A( \omega_k^j ) =\sum a_i \omega^{ij} = \sum a_i (w^{i\% k})^j = \sum a_i (\omega^j)^{i \%k} = B(\omega_k^j) \]

证毕。


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


*Bluestein

也叫 Chirp Z-Transform。

任意模长的 DFT，不要求 \(n\) 是 \(2^k\)。/qiang

曾经我们将卷积 \(\to DFT/IDFT\)，现在我们尝试把 \(DFT/IDFT\) 转化成卷积。

给定多项式 \(A(x) = a_0 + \cdots + a_{n-1} x^{n-1}\)，求其在 \(\omega^0,\cdots,\omega^{n-1}\) 的值，不要求 \(\omega\) 是单位根（可以是任意数，只要这是个等比数列就行）。

求 \(b_0 \sim b_{n-1}\)，其中 \(b_m = \sum_{k=0}^{n-1} a_k (\omega^m)^k\)。

多项式多点求值？太慢啦！


考虑重要的 trick

\[mk =\frac 1 2 \left [ m^{\underline 2} + (k+1)^{\underline 2} + (m-k)^{\underline 2} \right] \]

展开即可验证。

不能展开成 \(\frac 1 2 [m^k+k^2 - (m-k)^2]\)，在后面的推到中不保证 \(\omega\) 有二次剩余，所以非常局限！！！

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

\[b_m = \omega^{\frac 1 2 m^{\underline 2}} \sum_{i=0}^{n-1} a_k \omega^{\frac 1 2 (k+1)^{\underline 2} - \frac 1 2 (m-k)^{\underline 2}} \]

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


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

\[A_i = a_i \omega ^{\frac 1 2 (k+1)^{\underline 2}},B_{m+n} =b_m \omega^{\frac 1 2 m^{\underline 2}},C_j = \omega^{-\frac 1 2 (j-n)^{\underline 2}} \]

我们就有

\[B_{m+n} = \sum_{i=0}^{m+n} A_i C_{m+n-i} \]

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


对其进行一个推广。

求 \(A(x)\) 在 \(k\) 个 \(k\) 次单位根的值。

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


P6800 【模板】Chirp Z-Transform

板板。

给定一个 \(n\) 项多项式 \(P(x)\) 以及 \(c, m\)，请计算 \(P(c^0),P(c^1),\dots,P(c^{m-1})\)。所有答案都对 \(998244353\) 取模。

\(1\le n,m\le 10^6,0\le c,a_i<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 = \binom {m+k} 2 - \binom m 2 - \binom k 2 \]

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


*DFT 多路分治算法

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


我们先来考虑 \(n=3^k\) 时如何做 DFT？

给定 \(f(x) = a_0 + a_1 x + \cdots + a_{n-1} x^{n-1}\)，求它在 \(\omega^0,\cdots,\omega^{n-1}\) 的值。（\(\omega = \omega_n\)，是单位根）


发现我们可以类似 \(2^k\) 的做法，把整个式子分成三个部分

\[\begin{cases} a(x) = a_0 + a_3 x+ \cdots + a_{n-3} x^{\frac n 3 -1} \\ b(x) = a_1 + a_4 x + \cdots a_{n-2} x^{\frac n 3 -1} \\ c(x) = a_2 + a_5 x+ \cdots a_{n-1} x^{\frac n 3 -1} \end{cases} \]

我们定义 \(\sigma = \frac n 3\) 次单位根，那么有 \(\omega = \sigma ^{\frac 1 3},\sigma = \omega ^3\)。

于是只要我们先求出 \(a(x),b(x),c(x)\) 在 \(\sigma ^0,\cdots,\sigma^{\frac n 3 -1}\) 的值，就有

\[\begin{aligned} f(x) & = a(x^3) + xb(x^3) + x^2 c(x^3) \\ \therefore f(\omega^i) & = a(\omega ^{3i}) + x b(\omega^{3i}) + x^2 c(\omega^{3i}) \\ & = a(\sigma^i) + xb(\sigma^i) + x^2 c(\sigma^i) \end{aligned} \]

时间复杂度也是 \(\mathcal O(n \log n)\)。


同样地，我们可以把它推广到 NTT（也就是模 \(p\) 素数域）下面：

  • 卷积：\(n\) 放大到 \(3^k(n \le 3^k)\)，满足 \(3^k \mid p-1\)，取 \(\omega = g^{\frac {p-1} {3^k}}\)。
  • \(n\) 不能放大，\(n \mid p-1\)，且 \(n\) 的质因子较小，可以根据 \(n\) 的最小质因子分治下去。

*单位根反演

难道他真的是天才？

\[[ n \mid v] = \frac 1 n \sum_{j=0}^{n-1} (\omega_n^v)^j \]


证明是简单的：

由于

\[\sum_{j=0}^{n-1} x^j = \frac {1-x^n}{1-x} (x \neq 1) \]

设 \(x = \omega_n^v\)，则我们需要分两种情况讨论

  • \(x = 1\)，则 \(\frac 1 n \sum_{j=0}^{n-1} (\omega_n^v)^j=1\)，这种情况出现的条件是 \(n \mid v\)。
  • \(x \neq 1\)，则化简式的分母 \(1 - (\omega_n^v)^n = 1-(\omega_n^n)^v = 1-1=0\)，而分子不为 \(0\)，所以原式 \(=0\)，这种情况出现的条件是 \(n \nmid v\)。

得证。


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

\[\sum_v [n \mid v] F(v) \]

的式子。

P.S.：这里的 \(\omega\) 不一定是复数意义下的单位根，也可以是模 \(p\) 域的 \(n\) 次单位根。


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

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

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


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


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

\(T\) 组数据，每次给定 \(n,s,a_0,a_1,a_2,a_3\) 求

\[\sum_{i=0}^n \binom n i s^i a_{i \bmod 4} \mod 998244353 \]

\(1 \le T \le 10^5,1 \le n \le 10^{18},1 \le s,a_0,a_1,a_2,a_3 \le 10^8\)。

煎蛋题。


看到 \(\bmod\) 想到 单位根反演。

考虑答案为

\[\begin{aligned} & \sum_{t=0}^3 a_t \sum_{i=0}^n \binom n i s^i [4 \mid i-t] \\ = & \frac 1 4 \sum_{j=0}^3 \left( s \omega_4^j + 1 \right)^n \sum_{t=0}^3 a_t \omega_4^{-tj} \end{aligned} \]

结束！代码。

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)\times n\) 的有向图。这张图的每个顶点由一个二元组\((u,v)\)表示\((0\le u\le L,1\le v\le n)\)。
这张图不是简单图，对于任意两个顶点 \((u_1,v_1)(u_2,v_2)\)，如果 \(u_1<u_2\)，则从 \((u_1,v_1)\) 到 \((u_2,v_2)\) 一共有 \(w[v_1][v_2]\) 条不同的边，如果 \(u_1\ge u_2\) 则没有边。

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

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

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

问题来了：给定正整数 \(k\) 和 \(y\)（\(1\le y\le n\)），对于每个 \(t\)（\(0\le t<k\)），求有多少种舞曲（假设其长度为 \(m\)）满足 \(m \bmod k=t\)，且白兔最后停在了坐标第二维为 \(y\) 的顶点？

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

输出的结果对 \(p\) 取模。

\(p\) 为一个质数，\(10^8<p<2^{30}\)，\(1\le n\le 3\)，\(1\le x\le n\)，\(1\le y\le n\)，\(0\le w(i,j)<p\)，\(1\le k\le 65536\)，\(k\) 为 \(p-1\) 的约数，\(1\le L\le 10^8\)。


考虑走 \(m\) 步的方案数，容易发现是

\[\binom L m W^{m} 的 x 行y列的值 \]

于是，我们有

\[ans_t = \sum_{m=0}^L \binom L m W^m [k \mid m-t] \]

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

\[\begin{aligned} = & \frac 1 k \sum_{m=0}^L \binom L m W^m \sum_{j=0}^{k-1} \left( \omega_k^{m-t} \right)^j \\ = & \frac 1 k \sum_{j=0}^{k-1} \omega_k^{-tj} \sum_{m=0}^L \binom L m \left( \omega_k^j W \right)^m \\ = & \frac 1 k \sum_{j=0}^{k-1} \omega_k^{-tj} \left( I + \omega_k^jW \right)^L \end{aligned} \]

我们令 \(A_j = \left( I + \omega_k^jW \right)^L\)，\(a_j = A_j[x][y]\)，那么原式就是

\[\frac 1 k \sum_{j=0}^{k-1} \omega_{k}^{-tj} a_j \]

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


我们把 \(mk\) 展开成 \(\binom {m+k} 2 - \binom m 2 - \binom k 2\)，代码较简单，直接 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\),询问

\[\sum_{i=0}^n \binom n i \times p^{i} \times \left\lfloor \frac{i}{k} \right\rfloor \bmod 998244353 \]

\(1 \leq n,p <998244353,k \in \{2^{w}|0 \leq w \leq 20\}\)。

居然切掉了/jk


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

\[\begin{aligned} & \sum_{i=0}^n \binom n i \times p^i \times \left\lfloor \frac i k \right \rfloor \\ = & \sum_{i=0}^n \binom n i \times p^i \times \frac {i - i\bmod k} k \\ = & \frac 1 k \left[ \sum_{i=0}^n \binom n ii p^i - \sum_{i=0}^n \binom n i (i \bmod k ) p^i \right] \end{aligned} \]

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

\[\begin{aligned} & \sum_{i=0}^n \binom n i ip^i \\ = & \sum_{i=0}^n \binom {n-1}{i-1} n p^i \\ = & np \sum_{i=0}^{n-1} \binom {n-1} i p^i \\ = & np(p+1)^{n-1} \end{aligned} \]

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

\[\begin{aligned} & \sum_{i=0}^n \binom n i (i \bmod k) p^i \\ = & \sum_{t=0}^{k-1} t \sum_{i=0}^n \binom n i p^i [k \mid i-t] \\ = & \frac 1 k \sum_{t=0}^{k-1} t \sum_{i=0}^n \binom n i p^i \sum_{j=0}^{k-1} \left( \omega_k^{i-t} \right)^j \\ = & \frac 1 k \sum_{t=0}^{k-1} t \sum_{j=0}^{k-1} \omega_k^{-tj} \sum_{i=0}^n \binom n i \left(p \cdot \omega_k^j \right)^i \\ = & \frac 1 k \sum_{t=0}^{k-1} t \sum_{j=0}^{k-1} \omega_k^{-tj} \left(p \cdot \omega_k^j +1 \right)^n \end{aligned} \]

仔细一看，后面的 \(\sum_j\) 中不就是 \(\left < \left(p \cdot \omega_k^j +1 \right)^n \right >\) 的 IDFT 吗？

由于保证了 \(k = 2^{\omega}\)，所以可以直接 NTT。

于是就做完了，时间复杂度 \(\mathcal O(k \log k)\)。还跑得蛮快的。代码。

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,\cdots,n \}\)，她想知道有多少子集 \(T \subseteq S\)，满足 \(\gcd(\sum_{t \in T t,n}) = 1\)。

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

由于奈亚子设立了方便剧情发展的结界，所以你可以对 \(10^9+7\) 取模。

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

\(T \le 10,n \le 2 \times 10^{15}\)，保证 \(n\) 不为 \(10^9+7\) 的倍数。

乎✌的胜利！知乎。


令 \(m= \frac {n(n+1)} 2\)，那么问题即求

\[\sum_{i=1}^m [\gcd(i,n) = 1] f(i) \]

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

而我们发现 \(f(i)\) 的生成函数其实是 \(F(x) = \prod_{i=1}^n (1 +x^i)\)，也就有 \(f(i) = [x^i] F(x)\)。

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

\[\begin{aligned} & \sum_{i=1}^m [\gcd(i,n) = 1] f(i) \\ = & \sum_{i=1}^m \sum_{d \mid n,d \mid i} \mu(d) f(i) \\ = & \sum_{d \mid n} \mu(d) \sum_{i=1}^m [d \mid i] f(i) \end{aligned} \]

看到 \([d \mid i]\)，我们想打单位根反演，于是有

\[\begin{aligned} & \sum_{i=1}^m [d \mid i] f(i) \\ = & \frac 1 d \sum_{i=1}^m \sum_{j=0}^{d-1} \left ( \omega_d^i \right)^j f(i) \\ = & \frac 1 d \sum_{j=0}^{d-1} \sum_{i=1}^m \left( \omega_d^j \right)^i [x^i]F(x) \\ = & \frac 1 d \sum_{j=0}^{d-1} F \left( \omega_d^j \right) \end{aligned} \]

现在式子就变成了求

\[\sum_{d \mid n} \frac {\mu(d)} d \sum_{i=0}^{d-1} F(\omega_d^i) \]

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

\[\sum_{i=0}^{d-1} F(\omega_d^i) = \sum_{i=0}^{d-1} \left( \prod_{j=1}^n \left( 1+ \omega_d^{ij}\right) \right) \]

发现上面的 \(ij\) 是一个经典问题，也就是设 \(\gcd(i,d) = k\)，则 \(\omega_d^{ij}\) 会循环 \(k\) 次，每次循环节长度是 \(\frac d k\)，那么可以得到

\[\begin{aligned} & \sum_{i=0}^{d-1} \left( \prod_{j=1}^n \left( 1+\omega_d^{ij} \right) \right) \\ = & \sum_{k \mid d} \left( \sum_{i=0}^{d-1} [\gcd(i,d)=k] \right) \left( \prod_{i=1}^{ \frac d k } \left(1 + \omega_d^{ik} \right) \right)^{\frac {nk} d} \\ = & \sum_{k \mid d} \varphi \left( \frac d k \right) \left( \prod_{i=1}^{ \frac d k } \left(1 + \omega_d^{ik} \right) \right)^{\frac {nk} d} \end{aligned} \]

然后考虑结论

\[\prod_{i=1}^n (x - \omega_n^i) = x^n -1 \]

证明就是考虑 \(x^n-1\) 的 \(n\) 个根是 \(\omega_n^i\)，于是我们又可以将原式化简

\[\begin{aligned} & \sum_{k \mid d} \varphi \left( \frac d k \right) \left( \prod_{i=1}^{ \frac d k } \left(1 + \omega_d^{ik} \right) \right)^{\frac {nk} d} \\ = & \sum_{k \mid d} \varphi \left( \frac d k \right) \left( (-1)^{\frac d k} \prod_{i=1}^{ \frac d k } \left( -1 - \omega_d^{ik} \right) \right)^{\frac {nk} d} \\ = & \sum_{k \mid d} \varphi \left( \frac d k \right) \left( 1- (-1)^{\frac d k} \right)^{\frac {nk} d} \\ = & \sum_{k \mid d} \varphi \left( k \right) \left( 1 - (-1)^k \right)^{\frac n k} \\ = & \sum_{k \mid d,2 \nmid k} \varphi(k) 2^{\frac n k} \end{aligned} \]

于是答案就是

\[\begin{aligned} & \sum_{d \mid n} \frac {\mu(d)} d \sum_{k \mid d,2 \nmid k} \varphi(k) 2^{\frac n k} \\ = & \sum_{k \mid n,2 \nmid k} \varphi(k) 2^{\frac n k} \sum_{d\mid n,k \mid d} \frac {\mu(d)} d \end{aligned} \]

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

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

\[\begin{aligned} \mu(ij) & = \mu(i) \mu(j) [\gcd(i,j)=1] \\ \varphi(ij) & = \frac {\varphi(i) \varphi(j) \gcd(i,j)} {\varphi(\gcd(i,j))} \\ \to \varphi(i) & = \frac {\varphi(ij) \varphi(gcd(i,j))} {\varphi(j) \gcd(i,j)} \\ \sum_{d \mid n} \frac {\mu(d)} d & = \frac {\varphi(n)} n \end{aligned} \]

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

\[\begin{aligned} & \sum_{d \mid n,k \mid d} \frac {\mu(d)} d \\ = & \sum_{d \mid {\frac n k}} \frac {\mu(kd)} {kd} \\ = & \frac {\mu(k)} k \sum_{d \mid \frac n k} \frac {\mu(d)[\gcd(k,d)=1]} d \\ = & \frac {\mu(k)} k \sum_{d \mid \frac {n} {k \gcd(k,\frac n k)}} \frac {\mu(d)} d \\ = & \frac {\mu(k)} k \cdot \frac {k \gcd(k,\frac n k) \varphi \left(\frac n {k\gcd(k,\frac n k)} \right)} {n} \end{aligned} \]

设 \(x = k \gcd(k, \gcd(k,\frac n k))\)，于是

\[\begin{aligned} \varphi \left( \frac n x \right) = \frac {\varphi(n) \varphi(\gcd(\frac n x,x))} {\varphi(x) \gcd(\frac n x,x)} \end{aligned} \]

由于 \(\gcd(k\gcd(k,\frac n k),\frac {n} {k \gcd(k,\frac n k)}) = 1\)，所以

\[\varphi \left(\frac n x \right) = \frac {\varphi(n)} {\varphi(x)} \]

而我们还可以继续化简 \(\varphi(x)\)，也就是

\[\varphi(x) = \varphi \left(k \gcd\left(k, \frac n k \right) \right)=\frac {\varphi(k)\varphi \left ( \gcd \left(k,\frac n k \right) \right) \gcd\left(k, \gcd \left( k ,\frac n k \right) \right)} {\varphi \left(\gcd\left(k, \gcd \left( k ,\frac n k \right) \right) \right)} \]

稍微化简一下就是

\[\varphi(x) = \varphi(k) \gcd\left(k,\frac n k \right) \]

于是式子就可以继续化简

\[\begin{aligned} & \frac {\mu(k)} k \cdot \frac {k \gcd(k,\frac n k) \varphi \left(\frac n {k\gcd(k,\frac n k)} \right)} {n} \\ = & \frac {\mu(k)} k \cdot \frac {k\gcd \left(k,\frac n k \right) \varphi (n)} {n \varphi(k)\gcd\left(k,\frac n k \right)} \\ = & \frac {\varphi(n)} n \cdot \frac {\mu(k)} {\varphi(k)} \end{aligned} \]

那么

\[\begin{aligned} ans = & \sum_{k \mid n,2 \nmid k} \varphi(k) 2^{\frac n k} \sum_{d\mid n,k \mid d} \frac {\mu(d)} d \\ = & \frac {\varphi(n)} n \sum_{k \mid n,2 \nmid k} \varphi(k) 2^{\frac n k} \frac {\mu(k)} {\varphi(k)} \\ = & \frac {\varphi(n)} n \sum_{k \mid n,2 \nmid k} \mu(k) 2^{\frac n k} \end{aligned} \]

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

注意特判 \(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(x^a - 1, x^b - 1) + 1, x > 0\)。

特别的，如果 \(a = 0\) 或 \(b = 0\)，\(F(x, a, b) = 0\)。

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

令 \(L = F(m, a, b) + 1\)，\(R = F(m, c, d)\)。

问集合 \(\{L, L + 1, L + 2, \dots, R - 2, R - 1, R\}\) 有多少个子集和是 \(m\) 的倍数。

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

\(0 \lt m \le 10^7,L,R \le 10^{18},a,b,c,d \le 10^3,T \le 10^4\)。

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


考虑 \(L,R\) 究竟是什么？

对于 \(\gcd\)，我们不仅要考虑 辗转相除法，还要去思考 辗转相减法。

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

\[\gcd(x^{a} -1,x^{b}-1) = \gcd(x^{b}-1,x^{b} (x^{a-b}-1)) \]

然而，显然 \(x^b -1 \bot x^b\)，所以上式就变成了

\[\gcd(x^b-1,x^{a-b}-1) \]

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

\[\gcd(x^0-1,x^{p}-1) = x^{p}-1 \]

于是答案其实是

\[x^{\gcd(a,b)} \]

打表应该也能看出来。


那么现在我们知道了 \(L = m^{\gcd(a,b)}+1,R=m^{\gcd(c,d)}\)，容易得知 \(m \mid R-L+1\)，这一点会在之后的推到中有用。

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

首先列出答案式子

\[ans = \sum_{i=0}^{sum} [m \mid i] f(i) \]

其中 \(sum=\frac {(L+R)(R-L+1)} 2,f(i)\) 表示多少个集合的和是 \(i\)。

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

\[F(x) = \prod_{i=L}^R (1+x^i) \]

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

\[\begin{aligned} ans= & \sum_{i=0}^{sum} [m \mid i]f(i) \\ = & \frac 1 m \sum_{i=0}^{sum} \sum_{j=0}^{m-1} (\omega_m^{i})^j f(i) \\ = & \frac 1 m \sum_{i=0}^{m-1} F(\omega_m^i) \\ = & \frac 1 m \sum_{i=0}^{m-1} \left( \prod_{j=L}^R (1+\omega_m^{ij}) \right) \\ = & \frac 1 m \sum_{i=0}^{m-1} \left( \prod_{j=0}^{m-1} (1+ \omega_m^{ij}) \right)^{\frac {R-L+1} m} \\ = & \frac 1 m \sum_{d \mid m} \left( \sum_{i=0}^{m-1} [\gcd(i,m)=d] \right) \left( \prod_{j=0}^{\frac md} (1+\omega_{m}^{dj}) \right)^{d \cdot \frac {R-L+1} m} \\ = & \frac 1 m \sum_{d \mid m} \varphi(d) \left(1-(-1)^d \right)^{\frac {R-L+1} d} \\ = & \frac 1 m \sum_{d \mid m,2 \nmid d} \varphi(d) 2^{\frac {R-L+1} d} \end{aligned} \]

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


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

注意空集也要算，所以 \(L \gt 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(1 \le n \le 10^9,0 \le r \lt k \le 50,2 \le p \le 2^{30}-1)\)，求

\[\sum_{i \bmod k = r} \binom {nk} i \bmod p \]

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


直接变成 \((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;
}
posted @ 2024-07-23 15:46  H_W_Y  阅读(35)  评论(1)  编辑  收藏  举报