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。
以下的例题有些是没有放板子的,注意一下 清空 就可以了。但是每一遍都是重新写的。
有用的式子
证明两边同时求导即可。
下面例题中有些使用 vector
实现的多项式操作,具体存在于 CF755G 和 UOJ182。
*MTT
感觉还是需要单独提一嘴,后面单位根反演部分几乎是道道题要 MTT,真是令人难绷。/qd
其核心是把多项式系数
例如在 P4245 【模板】任意模数多项式乘法 中,乘出来的值高达
那么这道题就变成了求
直接做可以做到用
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 中的三步变两步,将其转化成复数多项式。
也就是我们可以设
那么就可以得到
于是我们需要的东西都可以从上面两个东西做加减出来了,容易发现这样只需要
果不其然会比上面快不少。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
给出
个数 ,定义 对
,求 的值。
, 。
典题先来一个。
首先推一下式子
我们设
于是我们有
但是这不是标准卷积形式,所以我们考虑 对称(镜像)。
令
这就是一个卷积形式了,直接 FFT 就可以了。
注意精度问题,我们需要给
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 - 容斥
的棋盘,有一些”超级车“,问:有多少个格子未被攻击。 一辆车可以攻击其所在行,列和对角线。
。
小容斥一手。
考虑把问题转化成求有多少个格子被至少一辆车攻击。
于是我们设
那么其实答案是
前面的
于是其实只用讨论如何求
不妨设
合理地对对角线编号,问题转化成求有多少
这种和式不就是 多项式卷积 形式?!
于是我们令
那么答案就是
这是可以直接 FFT 之后线性求出来的,结束!
时间复杂度
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_i+A_j+A_k = S i < j < k (i,j,k)$对数.
首先来考虑如果没有
容易发现其实就是
但是这里面有哪些重复的部分呢?
这个东西其实是
最后的部分是我们想要求的。
而对于第一部分,我们可以
再来分讨一下它是什么?
容易发现其实就是
这样通过两次 FFT 我们就得到了想要得到的东西了。
于是就做完了,时间复杂度
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 - 二维转一维
给定一张
的网格图。一些点会落在网格图上,记 为坐标是 的点的数量 ( )。现要求对于所有 成对的点之间的欧几里得距离,统计:
- 全部距离依次平方后,排序并去重后的结果
; - 上述序列去重前
重复出现的次数 ; - 所有距离的平均距离
。
, 。
应该是紫吧?
先写一个式子出来
这是二维的式子,非常不好处理,所以我们考虑一种映射关系
同时令
于是有
这个东西很像卷积,但是又不是标准的卷积形式。
按照之前的套路,我们直接对它 翻转,也就是令
于是就得到了
这就可以用卷积直接做了,上 FFT 就完事,不卡精度不卡空间(表示好评)。
注意特殊处理一下
时间复杂度
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
有一副超级扑克,包含无数张牌。对于每个正合数
,恰好有 张牌:黑桃 ,红桃 ,梅花 和方块 (分别用 、 、 和 表示)。没有其它类型的牌。 给定一个整数
,从 种花色中各选一张牌,问有多少种组合可以使得点数之和 ,例如, 时,有一种组合方法是 ,如下图所示。
不巧的是,有些牌已经丢失,并且为了让题目更毒瘤,我们还会提供两个正整数
和 ,你的任务是按顺序输出 时的答案。
。
唐唐板子题。
直接把四个多项式乘起来就可以了。
注意每组数据之后要输出一个 \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
有一排
个球,定义一个组可以只包含一个球或者包含两个相邻的球。现在一个球只能分到一个组中,求从这些球中取出 组的方案数。
, 。
首先,先写一个 dp,即设
于是我们有
由于需要
我的天,矩阵,快速,幂?
于是构造矩阵
那么答案其实就是
做完了,直接暴力 NTT 即可,常数大,时间复杂度
笔者本以为它过不了——但是它过了!需要开 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';
}
注意到矩阵乘法部分是可以优化的,我们发现可以先把
这样可以把常数减小一半。代码。
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。
考虑能不能从
发现是可以的
很明显的卷积形式,表示出来就是
于是我们尝试从
那么每次我们只需要维护
这就是 倍增 NTT 了,直接做还是 应该和上述代码类似),咕。
另外,这里还有一种
观察上面的式子
发现是一个递推方程的形式,所以我们考虑用 OGF 去优化它。
也就是设
于是根据递推形式可以得到
转化出来可以得到
我们需要求出
的形式。
于是直接列方程
带入
那么我们需要求的
注意到
最后的式子就可以以
这是一个 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 在字符串匹配中的应用
很久很久以前,在你刚刚学习字符串匹配的时候,有两个仅包含小写字母的字符串
和 ,其中 串长度为 , 串长度为 。可当你现在再次碰到这两个串时,这两个串已经老化了,每个串都有不同程度的残缺。 你想对这两个串重新进行匹配,其中
为模板串,那么现在问题来了,请回答,对于 的每一个位置 ,从这个位置开始连续 个字符形成的子串是否可能与 串完全匹配?
。
他都带通配符了,你就顺从它吧。
显然带着通配符,字符串的函数就不太能用了。
那怎么办呢?
考虑定义配对函数
当前字母是通配符时
考虑设
这并不好看。考虑按照之前的套路,我们令
这怎么这么像卷积形式?直接展开,可以得到
这三个部分都可以直接卷积处理了。
于是就做完了,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 在字符串匹配中的应用
给出一个门限值
和两个只包含 四种字符的基因串 和 。现在你要找出在下列规则中 在 中出现了几次。
在 的第 个位置中出现,当且仅当把 的首字符和 的第 个字符对齐后, 中的每一个字符能够在 中找到一个位置偏差不超过 的相同字符。 即对于所有的
,都存在一个 使得 且 。 例如
时, 出现在 的 号, 号和 号位置。 (编号从 开始。)
居然切了。
考虑按照上一题的思路,我们希望找到一个匹配函数来表示是否匹配这件事情。
由于 DNA 中只有四种碱基,所以我们考虑把它们每一个单独拿出来,构造一下两个序列
:若 则表示 位置合法,即存在一个和 距离不超过 的碱基;反之则 。 :若 则表示 是当前的碱基,反之则不是。
于是
就可以用 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 - 拉插
求
。
。
有些人拉插喜欢把
首先,为什么这是一个
一个序列做
因为每一次差分会让次数减少
由于我们把这个式子差分之后是
于是似乎就做完了,直接拉插
这东西稍微预处理一下就容易完成了。代码。
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
对于一个整数序列
,我们定义它的变换为 ,其中 ,其中 表示二进制按位或运算。 现在求有多少个长为
的序列 ,满足 ,使得它的变换 是严格单调递增的,对 取模。
, 。
发现多模 NTT 写成 vector 形式其实挺好写的。
首先,容易发现
写出 dp 方程,设
那么可以得到
看起来就非常卷积。
于是考虑倍增 NTT,考虑知道
就可以直接用多模 NTT 转移了。/qiang
时间复杂度
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
我的室友最近喜欢上了一个可爱的小女生。马上就要到她的生日了,他决定买一对情侣手环,一个留给自己,一个送给她。每个手环上各有
个装饰物,并且每个装饰物都有一定的亮度。 但是在她生日的前一天,我的室友突然发现他好像拿错了一个手环,而且已经没时间去更换它了!他只能使用一种特殊的方法,将其中一个手环中所有装饰物的亮度增加一个相同的非负整数
。并且由于这个手环是一个圆,可以以任意的角度旋转它,但是由于上面装饰物的方向是固定的,所以手环不能翻转。需要在经过亮度改造和旋转之后,使得两个手环的差异值最小。 在将两个手环旋转且装饰物对齐了之后,从对齐的某个位置开始逆时针方向对装饰物编号
,其中 为每个手环的装饰物个数, 第 个手环的 号位置装饰物亮度为 ,第 个手环的 号位置装饰物亮度为 ,两个手环之间的差异值为(参见输入输出样例和样例解释): 麻烦你帮他计算一下,进行调整(亮度改造和旋转),使得两个手环之间的差异值最小,这个最小值是多少呢?
, 。
居然自己做出来了 /jk
首先考虑把
根据之前字符串匹配的经验,和翻转的套路,这东西一看就非常卷积。
于是我们设
那么
简直就是卷积。于是可以 FFT 直接做。
而假设变化量是
于是就做完了,只需要一次 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
你需要求出
个点的简单 (无重边无自环) 有标号无向连通图数目。 由于这个数字可能非常大, 你只需要输出方案数对
( ) 取模即可。
。
有些人板子写挂调了一万年。
考虑经典结论,若
其中
于是设
时间复杂度
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
这个背包最多可以装
大小的东西 付公主有
种商品,她要准备出摊了 每种商品体积为
,都有无限件 给定
,对于 ,请你回答用这些商品恰好装 体积的方案数
, 。
考虑设
不知道怎么做?经典操作:先去 Ln,再去 Exp,可以得到
后面的式子是调和级数形式,只有
于是最后再做一次
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 - 多点求值
首先炮身的屏幕上显示了
个数,接着在模 一个质数意义下,给出了 条指令,每一条指令都是下列两种之一:
- 给所有数加上某一个数。
- 让所有数都变成原来的逆元。(保证这时所有数的逆元都存在)
在每个指令完成之后,你要回答当前屏幕上所有数的和。
跳蚤国王思索了片刻,发现这个问题奥妙重重,于是他让你——这附近最著名的民间科学家来帮他解决这个难题。
。
多项式写成 poly 还是要优秀一点/kk
我们设初始的序列为
而这里的
进而,我们可以将这个式子化简,根据初中数学的套路,该式可以被化简成
我们发现求和之后相当于是去求
的
而这个式子不是标准的多项式,看到它应该想到 求导!
对式子进行通分,那么分母就是
然而分子——是分母的 导数!这是可以用链式法则证明的。
于是就做完了,直接多点求值即可。代码。
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 有一棵树。这棵树上只有点
和 是黑色,其它的点都是白色。 每次,Tommy 可以将一个黑色
的点染成红色,然后把和 相邻的所有白色的点染成黑色。最后,所有的点都会被染成红色。 设第
个点是第 个被染成红色的,那么 是一个 到 的排列。Tommy 希望你帮他求出,有多少种不同的 。
。
一发过编过样例过题(怎么根 pkusc 场上一样),然后喜提最劣解(毕竟是取模)。
考虑不在
和 6.29 T3 一样,结论:随机排列满⾜树形拓扑序的概率是⼦树⼤小的倒数乘积。
证明可以去看那场联考的总结,写在里面了(7月训练记录)。
于是按照这个思路想下去,我们关心地其实就是
由于有两个点,我们尝试把它合成一个,也就是在建立一个虚点
那么现在就变成了基环树,我们只关心环上的点的顺序,而对于不在环上的点,它子树的贡献是一定的。
先给环上点标号,
枚举最后一个被染上的点?
怎么搞都是
由于是基环树——我们能不能考虑 断掉 一条边让他变成树?
发现是可以的,我们枚举基环树上的一条边
我们再钦定每个点的权值是
那么答案就是
其中
于是现在我们研究的其实是
想到 导数。发现这个式子其实是函数
在
根据 洛必达法则,我们对上下同时求导,也就是求
在
由于我们要同时求每个
于是就做完了,最终的答案式子就是
时间复杂度
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;
}
*降幂引理
根据名字就可以知道,这个引理用来解决指数
也就是我们现在有一个多项式
从
然而这两个多项式满足
即在
这是容易证明的,考虑直接把式子化开
证毕。
降幂引理 和 单位根反演 本质 相同!!!在之后的例题中我们会发现这一点。
*Bluestein
也叫 Chirp Z-Transform。
任意模长的 DFT,不要求
曾经我们将卷积
给定多项式
,求其在 的值,不要求 是单位根(可以是任意数,只要这是个等比数列就行)。 求
,其中 。
多项式多点求值?太慢啦!
考虑重要的 trick
展开即可验证。
不能展开成
,在后面的推到中不保证 有二次剩余,所以非常局限!!!
于是我们把这个式子带入之前的式子中
这不就是卷积形式?但是又不标准。
于是我们考虑一种非常简单的转化,令
我们就有
这就是一个标准卷积形式了,直接 FFT 或者 NTT 即可。常数可能有点大。
对其进行一个推广。
求
在 个 次单位根的值。
直接用降幂引理可以
P6800 【模板】Chirp Z-Transform
板板。
给定一个
项多项式 以及 ,请计算 。所有答案都对 取模。
.
没有预处理,直接暴力做!代码。
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;
}
似乎有更好的拆分
这可以避免一些神秘的负数问题。我也不知道是什么问题。
*DFT 多路分治算法
不如 Bluestein 应用广,但是不需要 FFT。
我们先来考虑
给定
,求它在 的值。( ,是单位根)
发现我们可以类似
我们定义
于是只要我们先求出
时间复杂度也是
同样地,我们可以把它推广到 NTT(也就是模
- 卷积:
放大到 ,满足 ,取 。 不能放大, ,且 的质因子较小,可以根据 的最小质因子分治下去。
*单位根反演
难道他真的是天才?
证明是简单的:
由于
设
,则 ,这种情况出现的条件是 。 ,则化简式的分母 ,而分子不为 ,所以原式 ,这种情况出现的条件是 。
得证。
单位根反演常常用来简化形如
的式子。
P.S.:这里的
有关 整除/同余 的题目 都可以 考虑 单位根反演!!!
有关 整除/同余 的题目 都可以 考虑 单位根反演!!!
有关 整除/同余 的题目 都可以 考虑 单位根反演!!!
上面扯了一堆看起来非常正确的东西,还是来练练题吧。
LOJ 6485. LJJ 学二项式定理 - 单位根反演
组数据,每次给定 求
。
煎蛋题。
看到
考虑答案为
结束!代码。
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] 白兔之舞 - 单位根反演
有一张顶点数为
的有向图。这张图的每个顶点由一个二元组 表示 。
这张图不是简单图,对于任意两个顶点,如果 ,则从 到 一共有 条不同的边,如果 则没有边。 白兔将在这张图上上演一支舞曲。白兔初始时位于该有向图的顶点
。 白兔将会跳若干步。每一步,白兔会从当前顶点沿任意一条出边跳到下一个顶点。白兔可以在任意时候停止跳舞(也可以没有跳就直接结束)。当到达第一维为
的顶点就不得不停止,因为该顶点没有出边。 假设白兔停止时,跳了
步,白兔会把这只舞曲给记录下来成为一个序列。序列的第 个元素为它第 步经过的边。 问题来了:给定正整数
和 ( ),对于每个 ( ),求有多少种舞曲(假设其长度为 )满足 ,且白兔最后停在了坐标第二维为 的顶点? 两支舞曲不同定义为它们的长度(
)不同或者存在某一步它们所走的边不同。 输出的结果对
取模。
为一个质数, , , , , , , 为 的约数, 。
考虑走
于是,我们有
直接单位根反演,再用二项式定理进行化简就可以得到
我们令
发现这是一个 IDFT 形式,于是可以直接用 Bluestein 做了。
我们把
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 小猪佩奇学数学 - 单位根反演
给定
,询问
。
居然切掉了/jk
很明显不能整除分块,于是我们考虑把下取整拆开
对于前者,我们考虑吸收公式
而对于后者,我们直接单位根反演
仔细一看,后面的
由于保证了
于是就做完了,时间复杂度
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)- 单位根反演
奈亚子很喜欢集合。
有一天真寻给了她一个集合
,她想知道有多少子集 ,满足 。 通俗地说,求集合
中,有多少子集的元素之和与 互素。 由于奈亚子设立了方便剧情发展的结界,所以你可以对
取模。 奈亚子因为可以生命时间加速,一瞬间就做出来了,而你只有
的时间,你还需要处理 组数据。
,保证 不为 的倍数。
乎✌的胜利!知乎。
令
其中
而我们发现
我们对这个式子进行莫比乌斯反演
看到
现在式子就变成了求
那么我们现在想去化简后面的部分,也就是
发现上面的
然后考虑结论
证明就是考虑
于是答案就是
我们考虑继续化简后面的那一块。
注意到一些数论函数的式子
都是可以用定义式或者计算式简单证明出来的,于是我们可以继续推导
设
由于
而我们还可以继续化简
稍微化简一下就是
于是式子就可以继续化简
那么
就做完了!实现非常简单。
注意特判
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 + 单位根反演
定义
。 特别的,如果
或 , 。 现在给出五个非负整数
。 令
, 。 问集合
有多少个子集和是 的倍数。 由于答案可能很大,你只需要输出方案数对
取模后的结果就可以了。
。
如果做了上一道题,那么这道题就跟坎瓜切菜一样。
考虑
对于
在这道题中,用 辗转相减 可以得到
然而,显然
发现系数也是更相减损!这也是一个 gcd 的过程,到最后
于是答案其实是
打表应该也能看出来。
那么现在我们知道了
于是类似于上一个题,我们就有了一下的推导。
首先列出答案式子
其中
容易知道
我们对原式做 单位根反演,可以得到
以上推导过程较为简略,但基本上和上一道题一样,看不懂的地方可以去看上一道题。
于是这样我们就得到了答案,直接枚举因数即可。
注意空集也要算,所以
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] 组合数问题 - 循环卷积
给定
,求
不知道为什么要用单位根反演推。
直接变成
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 中国大陆许可协议进行许可。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步