Poly 杂记（二）
由于放了代码，所以 md 的增长速度相当快。
为了让 cnblogs 的加载时间变短，就再开了一个。
书接上文。
CF438E The Child and Binary Tree - 开根 + 求逆
我们的小朋友很喜欢计算机科学，而且尤其喜欢二叉树。 考虑一个含有 \(n\) 个互异正整数的序列 \(c_1,c_2\cdots,c_n\)。如果一棵带点权的有根二叉树满足其所有顶点的权值都在集合 \(\{c_1,c_2,\cdots,c_n\}\) 中，我们的小朋友就会将其称作神犇的。
并且他认为，一棵带点权的树的权值，是其所有顶点权值的总和。
给出一个整数 \(m\)，你能对于任意的 \(1\leq s\leq m\) 计算出权值为 \(s\) 的神犇二叉树的个数吗？请参照样例以更好的理解什么样的两棵二叉树会被视为不同的。 我们只需要知道答案关于 \(998244353\) 取模后的值。
输入第一行有 \(2\) 个整数 \(n,m\) \((1\leq n,m\leq 10^5)\)。 第二行有 \(n\) 个用空格隔开的互异的整数 \(c_1,c_2\cdots,c_n\) \((1\le c_i\le10^5)\)。
输出 \(m\) 行，每行有一个整数。第 \(i\) 行应当含有权值恰为 \(i\) 的神犇二叉树的总数。请输出答案关于 \(998244353\) 取模的结果。
看到不连续的形如多项式的东西——考虑将它表示成连续，设 \(g_i\) 表示是否出现了 \(i\)。
再设 \(f_n\) 表示和为 \(n\) 的方案数，那么我们容易得到 dp 式子
其中 \(f_0 =1\)，这是一个标准的卷积形式啊，于是我们有
解出来
有两种情况，而取加号的时候 \(x \to 0\) 时该式 \(\to \infty\)，舍；取减号时则 \(\to 0\)，满足。
于是得到
但是 \(2G\) 显然没有逆——
于是化简一下，上下同时乘 \(1+\sqrt{1-4G}\) 得到
直接多项式开根 + 求逆即可。代码。
Code
poly Sqrt(poly F){
int I=qpow(2),P,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=A,B=Inv(B);
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*B[i]*I%H;
}
F=A,F.resize(n);
return F;
}
int main(){
/*2024.7.17 H_W_Y CF438E The Child and Binary Tree 多项式开根 + 求逆*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n>>m;G.resize(m+1);
for(int i=1,x;i<=n;i++){
cin>>x;
if(x<=m) G[x]=1;
}
for(int i=0;i<(int)G.size();i++) G[i]=(H-4ll*G[i]%H)%H;
G[0]=(G[0]+1)%H;
G=Sqrt(G);
G[0]=(G[0]+1)%H;
G=Inv(G);
for(int i=0;i<(int)G.size();i++) G[i]=2ll*G[i]%H;
for(int i=1;i<=m;i++) cout<<G[i]<<'\n';
return 0;
}
poly 版的求逆。
P3321 [SDOI2015] 序列统计 - 加与乘的转换
小C有一个集合 \(S\)，里面的元素都是小于 \(m\) 的非负整数。他用程序编写了一个数列生成器，可以生成一个长度为 \(n\) 的数列，数列中的每个数都属于集合 \(S\)。
小C用这个生成器生成了许多这样的数列。但是小C有一个问题需要你的帮助：给定整数 \(x\)，求所有可以生成出的，且满足数列中所有数的乘积 \(\bmod \ m\) 的值等于 \(x\) 的不同的数列的有多少个。
小C认为，两个数列 \(A\) 和 \(B\) 不同，当且仅当 \(\exists i \text{ s.t. } A_i \neq B_i\)。另外，小C认为这个问题的答案可能很大，因此他只需要你帮助他求出答案对 \(1004535809\) 取模的值就可以了。
\(1 \le n \le 10^9\)，\(3 \le m \le 8000\)，\(1\le x < m\)。
\(m\) 为质数，输入数据保证集合 \(S\) 中元素不重复。
什么操作可以完成加与乘之间的转换？指数函数和对数函数！
于是就做完了，直接找到 \(m\) 的原根即可。
注意由于我们需要循环卷积，所以写成暴力的快速幂形式。代码。
Code
poly operator *(poly F,poly G){
F.resize(P),G.resize(P);
NTT(F,P,1),NTT(G,P,1),Cro(F,G),NTT(F,P,-1);
for(int i=0;i<m-1;i++) F[i]=(F[i]+F[i+m-1])%H;
F.resize(m-1);return F;
}
poly Qpow(poly a,ll b){
poly res(m-1);res[0]=1;
while(b){if(b&1) res=res*a;a=a*a,b>>=1;}
return res;
}
int main(){
/*2024.7.17 H_W_Y P3321 [SDOI2015] 序列统计 加与乘的转换*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n>>m>>x>>len;
for(P=1;P<(m<<1);P<<=1);Init(P);
F.resize(m-1),init();
for(int i=2;i<m;i++) if(chk(i)){rt=i;break;}
pw[0]=1;
for(int i=1;i<m-1;i++) pw[i]=1ll*pw[i-1]*rt%m;
for(int i=0;i<m-1;i++) lg[pw[i]]=i;
for(int i=1,v;i<=len;i++){
cin>>v;v%=m;
if(!v) continue;
F[lg[v]]++;
}
F=Qpow(F,n);
cout<<F[lg[x]]<<'\n';
return 0;
}
CF1821F Timber - 多项式优化 dp
在 \([1, n]\) 的区间里放 \(m\) 棵树，每棵树的高度为 \(k\)。求有多少种放置树的方法，满足：
- 每个树都在整点上，且每个点最多只能放一棵树。
- 存在一种砍倒树的方案，使得树倒了之后不会越界，也不会有某个点被超过一棵树占据。你可以自由选择树向左倒（也就是占据区间 \([x - k, x]\)）或者向右倒（也就是占据区间 \([x, x + k]\)）。
\(1\leq n, m, k\leq 3\times 10 ^ 5\)。
不要跳步骤！可以从两个维度分别尝试！
考虑 dp，判断一种方案是否合法就是能往左倒往左倒，否则往右倒。
设 \(f_{i,j}\) 表示前 \(i\) 棵树，覆盖到 \(j\) 位置的方案数，枚举以下从哪里转移过来往左倒，哪里往右倒，于是可以得到
两者都是从 \(i-1\) 转移过来的，非常优秀！
于是我们考虑用多项式去刻画它们，设 \(F_i(x)\) 表示 \(\sum_{j=0}^{\infty} f_{i,j} x^j\)，于是我们有
由于 \(F_0(x)=1\)，所以我们有
答案就是
我们考虑继续化简 \(F_m(x)\)，容易发现，我们设
有
于是就做完了，\(O(n)\) 计算即可。代码。
Code
int main(){
/*2024.7.18 H_W_Y CF1821F Timber 多项式优化 dp*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n>>m>>k;init();
if(1ll*m*(k+1)>1ll*n) cout<<0<<'\n',exit(0);
for(int i=0;i<=m;i++){
int tmp=mul(mul(binom(m,i),pw[m-i]),binom(n-(m+i)*k,m));
if(i&1) del(ans,tmp);
else add(ans,tmp);
}
cout<<ans;
return 0;
}
ABC231G Balls in Boxes - EGF
有 \(n\) 个盒子，初始时第 \(i\) 个盒子内有 \(a_i\) 个小球，进行 \(k\) 次操作后，每次操作等概率随机选择一个盒子放入一个小球，设 \(k\) 次操作后每个盒子的小球个数为 \(b_i\)，那么得分为 \(\prod_{i=1}^n b_i\)。求出期望得分。
\(n\leq 1000,k,a_i\leq 10^9\)。
根本不用 FFT。
转化题目，其实是求
这种式子直接构造 EGF，设
那么我们要求的就是
稍微化简一下 \(F_i(x)\)，也就是把中间的加号拆开，那么
答案就是
前面可以直接展开，后面只需要暴力 \(O(n^2)\) 求。
\(m\) 很大，但只需要一个值，所以最后暴力卷积一下即可。代码。
Code
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define poly vector<int>
#define pb push_back
const int N=1e5+5,H=998244353;
int n,m,ans=0;
poly F,G;
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 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;
}
int main(){
/*2024.7.18 H_W_Y [ABC231G] Balls in Boxes EGF*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
F={1};cin>>n>>m;
for(int i=1,x;i<=n;i++){
G=F;cin>>x;
for(int i=0;i<(int)F.size();i++){
F[i]=mul(F[i],x);
if(i) add(F[i],G[i-1]);
}
F.pb(G.back());
}
for(int i=0,nw=1;i<min((int)F.size(),m+1);i++){
add(ans,mul(F[i],mul(qpow(n,m-i),nw)));
nw=mul(nw,m-i);
}
cout<<mul(ans,qpow(qpow(n,m)));
return 0;
}
AGC060D Same Descent Set - 容斥 + 多项式优化 dp
计算有多少个 \(1 \dots n\) 的排列对 \((A,B)\) 满足 \((A_{i+1}-A_i)(B_{i+1}-B_i) > 0\) 对于每一个 \(1\leq i < n\) 都成立。
\(2\leq n\leq 2\times 10 ^ 5\)。
好题啊。
对于一个排列，我们记录它 \(\lt\) 出现的位置集合 \(S \subseteq \{ 1,2,3,\cdots,n-1 \}\)，设 \(f(S)\) 表示满足这种关系的排列的个数，那么答案就是
考虑如何求 \(f(S)\)，发现这是一个经典问题，我们可以直接 容斥，设 \(g(S)\) 表示 \(S\) 中不等号未定义，其它都是 \(\gt\) 的方案数。
\(g(S)\) 是相对容易计算的，我们需要知道每一段的长度为 \(l_1 ,\cdots l_k\)，那么答案就是 \(\binom n {l1,l2,\cdots,l_k}\)。
于是这样我们可以容斥算出 \(f\)，也就是
那么答案就是
现在的问题来到 \(2^{|T_1\cup T_2|}\)，我们发现 \(|T_1 \cup T_2|\) 的自己数量是这个值，所以我们考虑直接用枚举自己来代替它，那么可以得到
现在回到最开始我们对 \(g\) 的计算，也就是一个集合 \(S\) 其实是把整个序列分成了 \(|S|+1\) 段极长的 \(\gt \cdots \gt\)。发现枚举 \(S\) 时已经把整个序列分成了若干段了。
那么我们现在考虑对这个东西 dp，过程中只需要记录上一次的端点在哪里就可以了。
也就是说，设 \(h_i\) 表示长度为 \(i\) 的序列分段方式的系数之和（最后一次在 \(i\)），就有
于是进而得到答案的 dp 方程，设 \(f_i\) 前 \(i\) 个元素 \(S\) 最后个元素在 \(i\) 的方案数，于是有
为什么没有 \(/ 4\)（当前这个元素放东西的 \((-2)^2\) 的贡献）？
我们发现 \(h\) 的计算过程中，我们是计算了最后一个元素的贡献的，所以在 \(f\) 的计算过程中是不用算的。
最后的答案就是 \(f_n \times 2^{n-1} \times (n!)^2 \times 4\)，这里 \((n!)^2\) 是 \(g\) 的计算中统一提出去的，而 \(4\) 是因为 \(n\) 位置本身不用算容斥系数，最后乘上就可以了。
观察一下，发现其实这两个递推式子都是卷积形式，几乎和分治 FFT 模板一模一样，于是稍微推一下既可以用多项式求逆/分治 FFT 做了。
时间复杂度 \(\mathcal O(n \log n)\) 或 \(\mathcal O(n \log^2n)\)。代码（用的多项式求逆）。
Code
poly calc(poly F){
for(int i=0;i<(int)F.size();i++) F[i]=(H-F[i])%H;
F[0]=(F[0]+1)%H;return Inv(F);
}
int main(){
/*2024.7.18 H_W_Y [AGC060D] Same Descent Set 多项式优化 dp + 容斥*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n;init();
F.resize(n+1);
for(int i=1;i<=n;i++) F[i]=mul(inv2,ifac[i]);
F=calc(F);
for(int i=1;i<=n;i++) F[i]=mul(F[i],F[i]);
F[0]=0,F=calc(F);
cout<<mul(mul(fac[n],fac[n]),mul(qpow(2,n+1),F[n]))<<'\n';
return 0;
}
CF553E Kyoya and Train - 分治 FFT
给定一张 \(n\) 个点 \(m\) 条边的无重边无自环的有向图，你要从 \(1\) 号点到 \(n\) 号点去。
如果你在 \(t\) 时刻之后到达 \(n\) 号点，你要交 \(x\) 元的罚款。
每条边从 \(a_i\) 到 \(b_i\)，走过它需要花费 \(c_i\) 元，多次走过同一条边需要多次花费。
走过每条边所需的时间是随机的，对于 \(k \in [1,t]\)，\(\frac{p_{i,k}}{10^5}\) 表示走过第 \(i\) 条边需要时间 \(k\) 的概率。因此如果多次走过同一条边，所需的时间也可能不同。
你希望花费尽可能少的钱（花费与罚款之和）到达 \(n\) 号点，因此每到达一个点，你可能会更改原有的计划。
求在最优决策下，你期望花费的钱数。
\(n \le 50\)，\(m \le 100\)，\(t \le 2 \times 10^4\)，\(x,c_i \le 10^6\)，\(\sum_{k=1}^t p_{i,k} = 10^5\)，答案精度误差 \(\le 10^{-6}\)。
蝴蝶变换数组开小了还报 TLE，真是勾勾又试试。
这和 THUSC 2022 的 D1 T1 非常像啊。
我们考虑倒着做 dp，设 \(f_{i,j}\) 表示到 \(i\) 的时间为 \(j\) 的期望，那么有
然后就可以转移了，时间复杂度 \(\mathcal O(mt^2)\)。
这时我们发现后面的那个式子非常像卷积，于是我们考虑化简它。
设 \(g_{i,j}\) 表示经过第 \(i\) 条边时间为 \(j\) 的期望，那么 \(f_{a_i,j} \leftarrow g_{i,j}\)，列出式子，发现
因为求 \(g\) 的时候需要后半部分的 \(f\)，所以我们考虑分治 FFT，考虑 \([mid+1,r]\) 对 \([l,mid]\) 的贡献。
我们发现 \(l \in [1,t],j+l \in [mid+1,r],j \in [l,mid]\)，于是转移式子其实是
直接 reverse 一下，设 \(a_j = f_{i,j+1},b_j=f_{b_i,r-j}\)，于是得到
就可以卷积了。
直接分治 FFT 即可，最开始还需要预处理 \(t \sim 2t\) 的值和两点之间的距离。代码。
Code
il void Mul(int n,int m){
int P;for(P=1;P<n+m;P<<=1);Init(P);
for(int i=n;i<P;i++) F[i]=0;
for(int i=m;i<P;i++) G[i]=0;
FFT(F,P,1),FFT(G,P,1);
for(int i=0;i<P;i++) F[i]=F[i]*G[i];
FFT(F,P,-1);
for(int i=0;i<n+m-1;i++) F[i]=F[i].real()/P;
}
il void solve(int l,int r){
if(l>=t) return;
if(l==r){
for(int i=1;i<n;i++) f[i][l]=inf;
for(int i=1;i<=m;i++) if(e[i].u!=n)
f[e[i].u][l]=min(f[e[i].u][l],g[i][l]+(ld)e[i].w);
return ;
}
int mid=((l+r)>>1);
solve(mid+1,r);
for(int i=1;i<=m;i++){
if(e[i].u==n) continue;
for(int j=0;j<r-l;j++) F[j]=p[i][j+1];
for(int j=0;j<r-mid;j++) G[j]=f[e[i].v][r-j];
Mul(r-l,r-mid);
for(int j=l;j<=mid;j++) g[i][j]+=F[r-j-1].real();
}
solve(l,mid);
}
int main(){
/*2024.7.18 H_W_Y CF553E Kyoya and Train 分治 FFT*/
scanf("%d%d%d%lf",&n,&m,&t,&val);
for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) if(i!=j) d[i][j]=inf;
for(int i=1;i<=m;i++){
cin>>e[i].u>>e[i].v>>e[i].w,d[e[i].u][e[i].v]=e[i].w;
for(int j=1;j<=t;j++) scanf("%lf",&p[i][j]),p[i][j]/=100000;
}
for(int k=1;k<=n;k++)
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
d[i][j]=min(d[i][j],d[i][k]+d[k][j]);
for(int i=0;i<2*t;i++) f[n][i]=i<=t?0:val;
for(int i=1;i<n;i++) for(int j=t;j<2*t;j++) f[i][j]=d[i][n]+val;
solve(0,2*t-1);
printf("%.10lf\n",f[1][0]);
return 0;
}
*自卷型分治 FFT
在上面的分治 FFT 中，我们都解决的是
自己卷别人 是相对容易的。
而现在我们来探究另一种分治 FFT：自己卷自己，也就是说
有什么区别呢？
我们发现原先的分治 FFT 中，每一次我们是把 \(f_{l \sim mid}\) 和 \(g_{1 \sim r-l+1}\) 卷起来，而这回，有可能 \(f_{1 \sim r-l+1}\) 没有准备好！
那怎么办呢？
考虑这种情况的出现 当且仅当 \(l=1\)，因为由于我们是分治，并且左区间一定 大于等于 右区间的长度，
所以如果 \(l \neq 1\)，那么 \(1 \sim r-l+1\) 一定是准备好了的。
那么唯一的问题就是 \(l=1\) 时——有些部分需要推迟。
推迟 如何处理，我们考虑每一次 \(l\neq 1\) 的左区间对右区间的贡献，这些贡献刚好是当时在某个 \(l=1\) 的区间中没有得到贡献的，
也就是那些推迟的贡献，但由于对称性，我们只需要把这部分贡献 \(\times 2\) 即可。
所以我在说什么？
其实自卷型分治 FFT 也就是两个部分：
- \(l=1\) 时，计算 \(f_{l \sim mid}\) 卷上 \(f_{l \sim mid}\)，贡献到 \(f_{mid+1,r}\)。
- \(l \neq 1\) 时，计算 \(f_{l \sim mid}\) 卷上 \(f_{1 \sim r-l+1}\)，将其 两倍 贡献到 \(f_{mid+1,r}\)。
于是就做完了，时间复杂度不变。\(\mathcal O(n \log^2n)\)。
UOJ 50 链式反应 - 分治 FFT
著名核物理专家 Picks 最近发现了一种非常适合核裂变的元素叫囧元素。囧元素原子序数为 1024，囧-2333 如果被一个中子撞击后会分裂成 蒟-1234 和 蒻-1098 同时释放出恰好 \(2\) 个中子，并且以破坏死光的方式释放光子。
核物理专家 Picks 做实验从来不用实验仪器。他用手指从宇宙中挑选出了 \(n\) 个 囧-2333 原子并编号为 \(1\) 到 \(n\)，并用帅气的眼神发射中子撞击了编号为 \(1\) 的囧原子，启动了链式反应。
当一个囧-2333原子被中子撞击时，有两种情况。要么这个囧-2333原子变为了囧-2334 并不再参与后续反应，要么囧-2333会进行核裂变，一方面，裂变产生的破坏死光会照射到 \(c\) 个不同的囧-2333原子并且这些原子会极为神奇地分解为氢和氦的混合物并不再参与后续反应。另一方面，裂变后产生的 \(2\) 个中子会分别撞上两个不同的囧-2333原子。显然被破坏死光照射后就不是囧-2333了，所以不可能又被中子撞上又被破坏死光照射到。
经过反复实验，核物理专家 Picks 终于确定了 \(c\) 的范围 \(A\)，其中 \(A\) 是一个非负整数集合。每次破坏死光会照射并影响到的囧-2333原子数量 \(c\) 满足 \(c \in A\) 。
链式反应时 Picks 会用肉眼观察实验中出现的事件（仅包括中子的撞击和破坏死光的信息），结束后 Picks 会写下实验记录。
但是很不幸 Picks 的实验记录丢失了。他只记得链式反应后没有剩余的囧-2333原子，而且每次一个囧-2333原子核裂变时，中子总是撞击编号比它大的囧-2333原子，破坏死光也总是照射编号比它大的囧-2333原子。
求可能会有多少种不同的实验记录。你需要对于 \(n = 1, \dots, n_{max}\) 输出答案。你只用输出答案对 \(998244353\)（\(7 \times 17 \times 2^{23} + 1\)，一个质数）取模后的值。
两个实验记录 \(T_1\)，\(T_2\) 被认为是相同的当且仅当：
- “编号为 \(v\) 的囧-2333分裂产生的中子撞上了编号为 \(u\) 的囧-2333” 这个事件在 \(T_1\) 中当且仅当这个事件在 \(T_2\) 中。
- “编号为 \(v\) 的囧-2333分裂产生的破坏死光照射了编号为 \(u\) 的囧-2333” 这个事件在 \(T_1\) 中当且仅当这个事件在 \(T_2\) 中。
\(n \leq 200000\)。
大清早上不去 UOJ 是什么啊/fn/fn/fn
自卷型分治 FFT。
考虑 dp，设 \(f_i\) 表示 \(i\) 点的答案，那么我们可以得到
由于两个子树无序，所以我们需要乘上 \(\frac 1 2\)。
发现式子涉及到 组合数，所以我们可以把它拆出来用 EGF 来完成，有
于是我们设 \(F_i = \frac {f_i} {i!}\)，那么有
容易相当定义 \(a_x = \frac {[x \in A]} {x!}\)，于是式子就变成了
三个数的卷积，如何处理？
考虑把两个东西合并一下，我们考虑合并 \(j,k\)，也就是枚举 \(j+k=s\)，设
那么
于是考虑 分治 FFT，同时计算 \(F,G\)。
设当前分治区间为 \([l,r]\)，目标是求 \(F_{l \sim r}\) 和 \(G_{l \sim r}\)：
- \(F_{l,mid}\) 去更新 \(G_{mid+1\sim r}\)，这是自卷型分治 FFT，有些操作需要 推迟。
- \(G_{l,mid}\) 去更新 \(F_{mid+1 \sim r}\)，直接做就可以了。
时间复杂度 \(\mathcal O(n \log^2 n)\)。
本质上根本不需要最后的转化，你可以直接三个一起卷。
否则常数巨大。不是编译原理大师，所以 vector 很慢不知道为什么。代码。
Code
#define mid ((l+r)>>1)
void solve(int l,int r){
// cerr<<"solve: "<<l<<','<<r<<endl;
if(l==r){
if(l==1) F[l]=1;
else F[l]=mul(F[l],qpow(2*l));
return;
}
solve(l,mid);
if(l==1){
int P;for(P=1;P<mid+1+mid+1+r;P<<=1);
Init(P);
for(int i=1;i<=mid;i++) a[i]=F[i];
NTT(a,P,1);
for(int i=0;i<P;i++) a[i]=a[i]*a[i]%H;
for(int i=0;i<r;i++) b[i]=A[i];
NTT(b,P,1);
for(int i=0;i<P;i++) a[i]=a[i]*b[i]%H;
NTT(a,P,-1);
for(int i=mid+1;i<=r;i++) (F[i]+=a[i-1])%=H;
}else{
int P;for(P=1;P<mid-l+1+r-l+r-l;P<<=1);
Init(P);
for(int i=l;i<=mid;i++) a[i-l]=F[i];
NTT(a,P,1);
for(int i=0;i<r-l;i++) b[i]=F[i];
NTT(b,P,1);
for(int i=0;i<r-l;i++) c[i]=A[i];
NTT(c,P,1);
for(int i=0;i<P;i++) a[i]=a[i]*b[i]%H*c[i]%H;
NTT(a,P,-1);
for(int i=mid+1;i<=r;i++) (F[i]+=mul(2,a[i-l-1]))%=H;
}
solve(mid+1,r);
}
void init(){
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=mul(fac[i-1],i);
ifac[n]=qpow(fac[n]);
for(int i=n;i>=1;i--) ifac[i-1]=mul(ifac[i],i);
}
int main(){
freopen("ex_reaction4.in","r",stdin);
freopen("1.out","w",stdout);
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n>>s;init();
for(int i=0;i<n;i++) if(s[i]=='1') A[i]=ifac[i];
solve(1,n);
for(int i=1;i<=n;i++) cout<<mul(F[i],fac[i])<<'\n';
return 0;
}
P4566 [CTSC2018] 青蕈领主 - 分治 FFT
小绿同学因为微积分这门课，对“连续”这一概念产生了浓厚的兴趣。小绿打算把连续的概念放到由整数构成的序列上，他定义一个长度为 \(m\) 的整数序列是连续的，当且仅当这个序列中的最大值与最小值的差，不超过\(m-1\)。例如 \(\{1,3,2\}\) 是连续的，而 \(\{1,3\}\) 不是连续的。
某天，小绿的顶头上司板老大，给了小绿 \(T\) 个长度为 \(n\) 的排列。小绿拿到之后十分欢喜，他求出了每个排列的每个区间是否是他所定义的“连续”的。然而，小绿觉得被别的“连续”区间包含住的“连续”区间不够优秀，于是对于每个排列的所有右端点相同的“连续”区间，他只记录下了长度最长的那个“连续”区间的长度。也就是说，对于板老大给他的每一个排列，他都只记录下了在这个排列中，对于每一个 \(1 \le i \le n\)，右端点为 \(i\) 的最长“连续”区间的长度 \(L_i\)。显然这个长度最少为 \(1\)，因为所有长度为 \(1\) 的整数序列都是连续的。
做完这一切后，小绿爬上绿色床，美美地做了一个绿色的梦。
可是第二天醒来之后，小绿惊讶的发现板老大给他的所有排列都不见了，只剩下他记录下来的 \(T\) 组信息。小绿知道自己在劫难逃，但是作为一个好奇的青年，他还是想知道：对于每一组信息，有多少个和信息符合的长度为 \(n\) 的排列。
由于小绿已经放弃治疗了，你只需要告诉他每一个答案对 \(998244353\) 取模的结果。
我们并不保证一定存在至少一个符合信息的排列，因为小绿也是人，他也有可能犯错。
\(1 \le T \le 100,1 \le N \le 50000,1 \le L_{i,j} \le j\)。
好题啊。
先来考虑如何判无解。
容易发现区间一定是 相离或包含，而最后一个数一定是 \(n\)。
因为只要区间相交，那么后面的区间一定就不优，和定义违背。
如果让每个区间向第一个 包含 它的区间连边，这样就会形成一棵树。
我们发现每个节点的贡献是一定的，也就是我们可以把它的每一个儿子看成一个点，那么这个长度为 \(cnt+1\) 的序列需要满足 任意一个连续的区间都包含最后一个节点。
我们不妨把这个记作 \(f_{cnt}\)（前面 \(+1\) 是因为还有它自己），那么其实答案就是
现在问题就变成了求 \(f_i\)。
考虑 dp，\(f_i\) 的定义是有多少个长度为 \(i+1\) 的序列满足 任意连续的区间都包含最后一个元素，也就是 \(L\) 的序列形如 \(\left<1,1,1, \cdots,1,n \right>\)。
由于每一次最后一个 元素 要变，所以我们考虑一种 映射，也就是将 位置 和 数值 翻转，这样对连续区间定义不变。
那么问题就变成了 任意一个连续的区间都包含最大值。
可是每一次 \(i+1\) 时又会改变最大值，这时不好维护的。
所以我们考虑在上一个 \(i\) 的基础上面插入 \(1\)，也就是让 \(i-1\) 的那些数从 \(2\) 开始，这样就非常好计数了。
我们分两种情况：
-
之前的序列 合法：
则当前插入的 \(1\) 只要不在 \(2\) 的两边就可以了，那么可以插的空隙数量是 \(i-1\)，贡献为 \((i-1)f_{i-1}\)。
-
之前的序列 不合法：
那么我们利用当前的这个 \(1\) 让序列合法，容易发现只会存在一个区间（可以不断包含）不合法。
因为如果存在两个 相离 的区间不合法，而只能插入一个 \(1\)，所以必然不能变为合法。
考虑枚举这个区间长度为 \(l\)，由于它不合法，所以一定不包含最大值；又因为插入 \(1\) 之后就合法了，所以一定不包含 \(2\)。
我们设这个区间的最小值是 \(x\)，那么区间值域是 \([x,x+l-1]\)，需要满足 \(x \gt 2,x+l-1 \le i,l \ge 2\)。
也就是说 \(l \in [2,i-2]\)，对于每一个 \(l\)，值域范围 \(x\) 有 \(i-l-1\) 中取值。
对于当前长度为 \(l\) 区间内的点，\(1\) 插入的合法位置满足 任意一个连续区间 都需要经过 \(1\)，我们发现这个和之前的 dp 定义非常像，也就是我们把 \(1\) 看成 \(l+1\)，那么其实合法的序列数量是 \(f_l\)。
进而，对于整体的贡献，我们就把这个长度为 \(l+1\) 的区间看成一个数，那么外面的答案就是 \(f_{i-l}\)。
所以这种情况的贡献就是 \(\sum_{l=2}^{i-2} (i-l-1) f_l f_{i-l}\)。
以上两种情况外面就得到了整个 dp 转移式
这很明显是一个 标准卷积形式，可以用 自卷型分治 FFT 解决。时间复杂度 \(\mathcal O(n \log^2n)\)。
注意到卷积中有 \((n-i-1)\) 的系数，于是每次对于 \(l \neq 1\) 的部分贡献时，我们的系数其实是 \((n-i-1)+(i-1)=n-2\)（之前推迟的贡献是 \(i-1\)），
于是可以直接卷积之后统一乘了。代码。
Code
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
const int N=1e6+5,H=998244353;
int n,m,U[N],p[N],tp=0,ans;
ull w[N],F[N],a[N],b[N];
struct Nod{
int l,r;
}st[N];
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(ull &a,int b){a=adc(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(int n){
for(int i=0;i<n;i++) U[i]=(U[i>>1]>>1)|((i&1)?n>>1:0),a[i]=b[i]=0;
}
void NTT(ull *p,int len,int op){}
#define mid ((l+r)>>1)
void solve(int l,int r){
if(l==r) return add(F[l],mul(F[l-1],l-1));
solve(l,mid);
if(l==1){
int P;for(P=1;P<mid+1+mid+1;P<<=1);Init(P);
for(int i=2;i<=mid;i++) a[i]=F[i];NTT(a,P,1);
for(int i=2;i<=mid;i++) b[i]=mul(F[i],i-1);NTT(b,P,1);
for(int i=0;i<P;i++) a[i]=a[i]*b[i]%H;
NTT(a,P,-1);
for(int i=mid+1;i<=r;i++) add(F[i],a[i]);
}else{
int P;for(P=1;P<mid-l+1+r-l+2;P<<=1);Init(P);
for(int i=l;i<=mid;i++) a[i-l]=F[i];NTT(a,P,1);
for(int i=2;i<=r-l+1;i++) b[i]=F[i];NTT(b,P,1);
for(int i=0;i<P;i++) a[i]=a[i]*b[i]%H;
NTT(a,P,-1);
for(int i=mid+1;i<=r;i++) add(F[i],mul(i-2,a[i-l]));
}
solve(mid+1,r);
}
void SOLVE(){
ans=1,tp=0;
bool fl=true;
for(int i=1,ct,x;i<=n;i++){
ct=0;cin>>x;
int l=i-x+1,r=i;
while(tp&&st[tp].r>=l){
if(st[tp].l<l) fl=false;
++ct,--tp;
}
st[++tp]={l,r},ans=mul(ans,F[ct]);
}
if(!fl||tp!=1) return cout<<0<<'\n',void();
cout<<ans<<'\n';
}
int main(){
/*2024.7.22 H_W_Y P4566 [CTSC2018] 青蕈领主 分治 FFT*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
int _;cin>>_>>n;
F[1]=2,F[0]=1,solve(1,n);
while(_--) SOLVE();
return 0;
}
*常系数齐次线性递推
此部分内容涉及线性代数相关知识，但目前这里并没有太多证明。
一个序列 \(f_i\) 若满足
\[f_i = c_1 f_{i-1}+c_2 f_{i-2}+\cdots +c_k f_{i-k} \]其中 \(c_i\) 是常数，我们称其为 \(k\) 次常系数递推式。
给定 \(f_{0 \sim k},c_{0 \sim k},k \le 32000,n \le 10^9\)，如何快速求出 \(f_n\) 呢？
于是我们引入 FFT 优化常系数齐次递推式。
首先，看到这个式子，我们容易想到 矩阵快速幂，也就是构造出矩阵
那么
于是直接做时间复杂度 \(\mathcal O(k^3 \log n)\)。
而我们的目标是去求 \(M^n\)。
如何优化，我们引入 \(M\) 的 特征多项式 为
这是可以用 \(\det(\lambda I - M)\) 的 \(A_{1,*}\) 计算出来的。
根据 Cayley - Hamilton 定理，我们有
于是只要我们构造出一个 \(k-1\) 次多项式 \(r(x)\) 满足
也就是多项式取模操作。
如果我们把 \(x=M\) 带进去就可以得到
设
那么
我们只需要 \(f_n\) 的值，也就是说
只要我们求出了 \(r(x)\)，那么 \(f_n\) 就可以计算了。
现在问题就变成了求
由于 \(n\) 太大了，但是考虑 \(\bmod\) 的性质：
于是就可以用多项式取模倍增做啦，时间复杂度 \(\mathcal O(k \log k\log n)\)。
可以通过 P4723 【模板】常系数齐次线性递推！代码。
Code
poly operator %(poly F,poly G){
int n=(int)F.size(),m=(int)G.size(),P;
if(n<m) return F;
poly A,B;P=n-m+1;
A=F,Rev(A),A.resize(P);
B=G,Rev(B),B.resize(P);
Inv(B),B=A*B,B.resize(P),Rev(B),B=B*G;
for(int i=0;i<m-1;i++) G[i]=(F[i]-B[i]+H)%H;
G.pop_back();return G;
}
poly Qpow(int n,poly mod){
poly res={1},a={0,1};
res=res%mod,a=a%mod;
while(n){if(n&1) res=(res*a)%mod;a=(a*a)%mod,n>>=1;}
return res;
}
int main(){
/*2024.7.20 H_W_Y P4723 【模板】常系数齐次线性递推*/
cin>>n>>k;
F.resize(k+1);F[k]=1;
for(int i=k-1,x;i>=0;i--){
cin>>x;x=(x%H+H)%H;
F[i]=(H-x)%H;
}
F=Qpow(n,F);
for(int i=0,x;i<k;i++){
cin>>x;x=(x%H+H)%H;
add(ans,mul(x,F[i]));
}
cout<<ans;
return 0;
}
*多项式平移
多项式 \(f(x) = \sum_{i=0}^n a_i x^i\)，求 \(f(x+c)\)。
多项式平移本质上属于多项式的复合，但这是最简单的版本。
我们直接将式子展开，可以得到
看起来就非常卷积，于是我们设
也就可以得到
发现这是非标准形式，于是我们考虑 reverse 一下，令
于是有
这就是 标准的卷积形式 了，直接做就可以做到 \(\mathcal O(n \log n)\)。
*下降幂多项式（FFP）
一个 \(n\) 次多项式可以唯一表示成 \(f(x) = \sum_{i=0}^n b_i x^{\underline i}\)。
这样的多项式我们就称之为下降幂多项式。
求和/斯特林数 问题都可以考虑 下降幂多项式 解决。
接下来我们来探讨一些下降幂多项式的操作。
已知 \(f(x),\left<b_i \right>\)，求 \(f(x+c)\) 的下降幂表示。
首先，我们知道下降幂也有 二项式定理，也就是
于是类似于上面普通幂多项式平移的方法，我们直接展开就可以得到
同样还是设
就可以卷积了。和上面几乎一模一样，时间复杂度 \(\mathcal O(n \log n)\)。
而下降幂多项式，有一个非常好的性质的就是
*多项式连续点值 \(\longleftarrow {\mathcal O(n \log n)} \longrightarrow\) 下降幂表示的系数
这也就是它广泛应用的一大原因。
如何完成这个转化，我们先考虑一种特殊情况。
假定 \(f(x)\) 的 \(0 \sim n\) 的点值为 \(y_0 \sim y_n\)，我们需要它和 \(f(x) = \sum_{i=0}^n b_i x^{\underline i}\) 互相转化。
注意：当 \(k \lt i\) 时，\(k^{\underline i} =0\)。
我们分两种情况来完成：
-
已知后者推前者：
把 \(y_k\) 表示出来，也就是
\[y_k = f(k) = \sum_{i=0}^n b_i k^{\underline i} =\sum_{i=0}^k b_i \frac {k!} {(k-i)!} \]那么
\[\frac {y_k} {k!} = \sum_{i=0}^k b_i \frac{1} {(k-i)!} \]发现这是一个 标准卷积形式，直接卷积即可完成。
-
已知后者推前者：
我们令 \(\hat y(x) = y_0 + \frac {y_1} {1!} x + \frac {y_2} {2!} x^2 + \cdots\)，那么上式就说明了
\[\hat y (x) = b(x) e^{x} \pmod {x^{n+1}} \]那么
\[b(x) = \hat y(x) e^{-x} \pmod {x^{n+1}} \]于是又可以卷积了。
时间复杂度都是 \(\mathcal O(n \log n)\)。
现在我们考虑其它的情况。
\(f(x)\) 在 \(c \sim n+c\) 的点值 \(y_0 \sim y_n\)，和 \(\left< b_i \right>\) 间的相互转化。
考虑直接构造多项式平移 \(g(x) = f(x+c)\)，
那么问题就变成了求 \(g(x)\) 的 \(0 \sim n\) 的点值，于是就可以 \(\mathcal O(n \log n)\) 做了。
总结一下：
我们发现普通幂转下降幂需要 \(\mathcal O(n \log ^2n)\) 的时间，瓶颈在于 多点求值，求普通幂连续点的点值。
而接下来从连续点值转化到下降幂就只需要 \(\mathcal O(n \log n)\) 的时间了。
而 普通幂，做 等比数列 的求值可以做到 \(\mathcal O(n \log n)\)；
下降幂，做 连续点值 可以做到 \(\mathcal O(n \log n)\)。
Trick：连续点值求和 用 下降幂！！！
P5394 【模板】下降幂多项式乘法
给定一个 \(n\) 次下降幂多项式 \(A(x)\) 和 \(m\) 次下降幂多项式 \(B(x)\)，你要求出一个 \(n+m\) 次下降幂多项式 \(F(x)\) 满足 \(F(x)=A(x)B(x)\)。
由于结果会很大，你输出的多项式的系数应对 \(998244353\) 取模。
\(1\leqslant n,m\leqslant 10^5\)，\(a_i,b_i\in[0,998244353)\)，\(a_n,b_m \neq 0\)。
直接用上面的 trick，我们先求出 \(0 \sim n+m\) 的点值，在转回去即可。
时间复杂度 \(\mathcal O(n \log n)\)。代码。
Code
poly FDT(poly F,int op){
int n=(int)F.size();
poly a;a.resize(n);
if(op==1) for(int i=0;i<n;i++) a[i]=ifac[i];
else{
for(int i=0;i<n;i++){
if(i&1) a[i]=dec(0,ifac[i]);
else a[i]=ifac[i];
F[i]=mul(F[i],ifac[i]);
}
}
F=F*a;F.resize(n);
if(op==1) for(int i=0;i<n;i++) F[i]=mul(F[i],fac[i]);
return F;
}
int main(){
/*2024.7.22 H_W_Y P5394 【模板】下降幂多项式乘法*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n>>m;++n,++m;
F.resize(n+m-1),G.resize(n+m-1);
for(int i=0;i<n;i++) cin>>F[i];
for(int i=0;i<m;i++) cin>>G[i];
init(n+m);
F=FDT(F,1),G=FDT(G,1);
for(int i=0;i<n+m-1;i++) F[i]=F[i]*G[i]%H;
F=FDT(F,-1);
for(int i=0;i<n+m-1;i++) cout<<F[i]<<' ';
return 0;
}
CF865G Flowers and Chocolate - 常系数线性递推
给定长为 \(F\) 的序列 \(p\)、长为 \(B\) 的序列 \(c\) 和整数 \(N\)，你需要对满足如下条件的二元组 \((f,b)\) 计数：
- \(f\) 是由 \(1\) 至 \(F\) 的整数组成的长为 \(N\) 的序列，\(b\) 是由 \(1\) 至 \(B\) 的整数组成的序列；
- \(\sum\limits_{i=1}^Np_{f_i}=\sum\limits_{i=1}^{|b|}c_{b_i}\)，其中 \(|b|\) 代表序列 \(b\) 的长度。
答案对 \(10^9+7\) 取模。
牛牛题，好极了。
首先考虑生成函数，容易发现
发现其实我们可以把 \(G(x)\) 表示成 常系数齐次线性递推形式，也就是设 \(t_i\) 表示 \(i\) 出现的次数，那么有特征多项式
其中 \(m\) 是最大值。
但是还是发现 \([x^i] G(x)\) 不好求，于是我们考虑 把 \(G(x)\) 向右平移 \(m-1\) 位，这样前 \(m\) 位只有 最后一位是 \(1\) 了！！！
所以
那么答案就是
于是就做完了。
并不用 MTT，可以直接 \(n^2\) 做多项式乘法和取模，取模时则是抓住 \(p[m]=1\) 的性质，每次减少一次，这样不断算下去的。
这样我们就做完了，时间复杂度 \(\mathcal O(m^2 \log n)\)。代码。
Code
#include <bits/stdc++.h>
using namespace std;
#define poly vector<int>
#define ll long long
const int N=1e3+5,H=1e9+7;
int n,m,ct[N],k,p[N];
ll len;
poly F,P,G;
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);}
poly operator *(poly a,poly b){
poly r(m<<1);
for(int i=0;i<m;i++) for(int j=0;j<m;j++) add(r[i+j],mul(a[i],b[j]));
for(int i=(m<<1)-1;i>=m;i--) for(int j=m-1;j>=0;j--) del(r[i+j-m],mul(r[i],P[j]));
r.resize(m);return r;
}
poly operator +(poly a,poly b){
for(int i=0;i<m;i++) add(a[i],b[i]);
return a;
}
poly qpow(poly a,ll b){
poly res(m);res[0]=1;
while(b){if(b&1) res=res*a;a=a*a,b>>=1;}
return res;
}
int main(){
/*2024.7.26 H_W_Y CF865G Flowers and Chocolate 常系数齐次线性递推方程*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n>>k>>len;
for(int i=1;i<=n;i++) cin>>p[i];
for(int i=1,x;i<=k;i++) cin>>x,m=max(m,x),++ct[x];
m=max(m,2);
P.resize(m+1);P[m]=1;
for(int i=1;i<=m;i++) del(P[m-i],ct[i]);
F.resize(m),G.resize(m),G[1]=1;
for(int i=1;i<=n;i++) F=F+qpow(G,p[i]);
F=qpow(F,len);
G[1]=0,G[m-1]=1,F=G*F;
cout<<F[m-1]<<'\n';
return 0;
}
CF1096G Lucky Tickets - GF
一个 \(n\) 位数，每位可以是给出的 \(k\) 个数码中的一个数，可以有前导 \(0\)，输出前 \(\frac n 2\) 位之和与后 \(\frac n 2\) 位之和相等的方案数，保证 \(n\) 是偶数。
输入的第一行是两个整数 \(n,k\)，接下来的一行有 \(k\) 个数 \(d_1,d_2,\cdots,d_k(0\leq d_i\leq 9)\)。
输出一个数，为方案数模 \(998244353\) 的值。
\(2 \le n \le 2 \times 10^5,k \le 10\)。
tang。
你生成函数一下，就变成了求
于是可以直接算 点值，直接做完了。
甚至还可以又简单做法
这样可以 \(\mathcal O(nk^2)\) 去递推了。
但是我写的是暴力做法，\(\mathcal O(nk \log n)\)。代码。
Code
int main(){
/*2024.7.26 H_W_Y CF1096G Lucky Tickets tang*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n>>m;n/=2;
F.resize(10);
for(int i=1,x;i<=m;i++) cin>>x,F[x]++;
int P;for(P=1;P<n*10;P<<=1);Init(P);
F.resize(P),NTT(F,P,1);
for(int i=0;i<P;i++) F[i]=qpow(F[i],n);
NTT(F,P,-1);
for(auto i:F) add(ans,mul(i,i));
cout<<ans;
return 0;
}
*IDFT 拓展
给定 \(q^0,\cdots,q^{n-1}\) 各不相等，\(f\) 为 \(n\) 次多项式，\(y_i = f(q^i)\)。
在 Poly 杂记（一）中，我们用 bluestein 算法实现了给定多项式 \(f\)，对每个 \(y_i\) 的求值。
而在这里，我们考虑直接变成 IDFT 形式，及给定 \(y\)，求 \(f\)。
乍一看似乎可以用Bluestein 继续解决，但实际上，我们推一推会发现 Bluestein 需要用到 \(y_0 \sim y_{2n-1}\) 的值，后面的 \(n\) 个值我们是不知道的，非常遗憾用不了。那么现在我们来讲另外一种方法，过程中会用到 Bluestein，但是并不是那么直接。
和快速插值一样，我们还是先考虑 拉格朗日插值，设
那么根据拉格朗日插值，我们有
这是用拉插最简单的式子进行一下 洛必达法则 就可以得到的。
那么现在步骤分为一下三步：
-
算 \(A(x)\)：考虑用分治 FFT 解决。
以 \(A(x) = (x-q^0) (x-q^1) (x-q^2) (x-q^3)\) 为例，假设我们算出了
\[h(x) = (x-q^0)(x-q^1) \]那其实对于后半部分，是可以直接线性推出来的，也就是
\[\begin{aligned} h\left(\frac x {q^2} \right) = & \left( \frac {x} {q^2} - q^0 \right) \left(\frac x {q^2} - q^1 \right) \\ = & \frac 1 {q^2} \frac 1 {q^2} (x-q^2) (x-q^3) \end{aligned} \]这就可以直接递推了。
所以这样的分治 FFT 的时间是
\[T(n) = T(n/2) + n\log n = \mathcal O(n \log n) \] -
算 \(A'(x)\) 在 \(q^0,q^1,\cdots,q^{n-1}\) 的值，这是一个标准的 Bluestein 形式，可以直接完成，时间复杂度 \(\mathcal O(n \log n)\)。
我们记 \(d_i = \frac {y_i} {A'(q^i)}\)。
-
算 \(\sum_{i=0}^{n-1} \frac {d_i} {x-q^i}\)。
考虑对 \(\frac 1 {x-q^i}\) 进行一些简化，以至于它可以被展开，也就可以得到
\[\frac 1 {x-q^i} = - \frac {q^{-i}} {1-q^{-i} x} \]于是这个式子其实就是
\[\sum_{i=0}^{n-1} \sum_{j=0} q^{-i(j+1)} x^j d^i = \sum_{j=0}^{n-1} \left( \sum_{i=0}^{n-1} q^{-i(j+1)} d_i \right)x^j \]如果我们定义
\[D(x)= \sum_{i=0}^{n-1} d_i x^i \]那么我们求的其实就是 \(D(x)\) 在 \(q^{-1},q^{-2},\cdots,q^{-n}\) 的值，这也是可以用 Bluestein 直接解决的。时间复杂度也是 \(\mathcal O(n \log n)\)。
上述三个步骤的时间复杂度都是 \(\mathcal O(n \log n)\)，所以总时间复杂度也是 \(\mathcal O(n \log n)\) 的。
看起来常数就挺大的。
*数幂求和
非常牛牛的一部分。
Part 1
给定 \(n,m\)，求 \(s_0,\cdots,s_m\)，其中 \(s_i = 0^i + 1^i + \cdots + (n-1)^i\)。
这是 数幂求和 的简单版本，我们之后再给出一个更为复杂的。
容易发现这是可以用 EGF 直接求解的，设
即 \(s\) 数组的 EGF。那么我们把 \(s_i\) 展开就可以得到
最后一步用到了等比数列求和，这似乎就可以用多项式求逆解决了。但是我们注意到首项并不为 \(1\)，所以其实我们还可以化简。
也就是分子是
同理，分母是
于是答案就是
直接多项式求逆的时间复杂度是 \(\mathcal O(m \log m)\)。
而这个东西其实和 伯努利数 有一定的联系，笔者还没有去研究这个东西，也就是说由于我们知道伯努利数的 EGF 是
然而这里的 EGF 是
所以这里是有一些联系的，进而，我们有如下的公式
证明先咕了。
其实我们还有另外一种方法，首先，将 \(n\) 看成 \((1+(n-1))\)，于是可以得到
把这些式子都加起来，那么左边抵消掉了。于是就有
然后两边同时 \(\div (i+1)!\) 就可以得到
这很明显是一个卷积形式，于是这个狮子相当于
也就和上面一种方法推出来一样了。
时间复杂度 \(\mathcal O(n \log n)\)。
Part 2
我们对之前的问题进行一个加强
给定 \(n,m\)，\(a_0,a_1,\cdots,a_{n-1}\)，求 \(s_0,\cdots,s_m\)，其中 \(s_i = {a_0}^i + {a_1}^i + \cdots {a_{n-1}}^i\)
用上面的 EGF 的方法，我们发现最终画出来的式子并不是一个 等比数列 形式，所以是不太可以做的。
于是我们考虑 OGF，定义 \(S(x)\) 为 \(s_0 \sim s_m\) 的 OGF，那么我们有
最后的式子怎么算？我们考虑用 分治 FFT 来完成。
定义
那么我们需要维护的其实是
那么每一层转移时我们有
于是时间复杂度是 \(\mathcal O(n \log^2n)\) 的。
Part 3 - 经典应用
给定多项式 \(H(x) = \sum_{j=0}^{m-1} h_j x^j\)，计算 \(H(a_0x) + H(a_1 x) + \cdots + H(a_{n-1}x)\)。
还是比较简单吧。
直接展开式子
对于每个 \(x^j\) 的式子，其实是我们在 Part 2 中的问题，于是可以 \(\mathcal O(n \log^2 n)\) 解决。
BZOJ3284 不等式
UOJ 310 黎明前的巧克力 - FWT
有 \(n\) 个巧克力，每个有属性 \(a_i\)。你需要从中选出两个不交的集合，使得这两个集合权值相同（权值为集合内巧克力属性的异或和）。
问方案数对 \(998244353\) 取模的余数。
\(1 \le n \le 10^6,1 \le a_i \le 10^6\)。
nit orz。
首先我们容易发现一个 dp，即设 \(dp_{i,j}\) 表示前 \(i\) 个巧克力，异或和为 \(j\) 的方案数。
那么容易得到转移方程
于是答案其实是求
其中 \(\prod\) 是 异或卷积。于是我们就有了 \(n\) 次 FWT 的做法。
显然不优，但是你考虑对这个东西进行 FWT，每次出来的这一位一定要么是 \(3\) 要么是 \(-1\)。
由于我们知道 FWT 的线性性：多个 FWT 的和等于和的 FWT。
于是考虑把 \((1+2x^{a_i})\) 都加到一起然后一起 FWT，于是每一位形如 \(3a-b\) 其中 \(a,b\) 还满足 \(a+b=n\)。
稍微解一下方程就可以得到最终的答案，再 FWT 回去就可以了。
于是我们就做完了，时间复杂度 \(\mathcal O(n \log n)\)。代码。
Code
#include <bits/stdc++.h>
using namespace std;
const int N=2e6+5,H=998244353;
int n,m,inv2,a[N],lim=1;
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 FWT(int *p,int lim,int op){
for(int le=1;le<lim;le<<=1)
for(int i=0;i<lim;i+=(le<<1))
for(int j=i;j<i+le;j++){
int x=p[j],y=p[j+le];
if(op) p[j]=adc(x,y),p[j+le]=dec(x,y);
else p[j]=mul(adc(x,y),inv2),p[j+le]=mul(dec(x,y),inv2);
}
}
int main(){
/*2024.8.19 H_W_Y UOJ310. 【UNR #2】黎明前的巧克力 集合幂级数*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n;inv2=qpow(2);
for(int i=1,x;i<=n;i++){
cin>>x,a[0]++,a[x]+=2;
while(lim<=x) lim<<=1;
}
FWT(a,lim,1);
for(int i=0;i<lim;i++){
int A=mul(adc(a[i],n),qpow(4)),B=mul(dec(3*n,a[i]),qpow(4));
a[i]=mul(qpow(3,A),qpow(H-1,B));
}
FWT(a,lim,0);
cout<<dec(a[0],1);
return 0;
}
ABC367G Sum of (XOR^K or 0) - FWT
给你正整数 \(N\) ， \(M\) ， \(K\) 和一个非负整数序列： \(A=(A_1,A_2,\ldots,A_N)\) 。
对于一个非空的非负整数序列 \(B=(B_1,B_2,\ldots,B_{|B|})\) ，我们对其分数定义如下。
- 如果 \(B\) 的长度是 \(M\) 的倍数： \((B_1 \oplus B_2 \oplus \dots \oplus B_{|B|})^K\)
- 否则 \(0\)
这里， \(\oplus\) 表示位向 XOR。
求以 \(998244353\) 为模数， \(A\) 的 \(2^N-1\) 个非空子序列的分数之和。
什么是位向 XOR？非负整数 \(A\) 和 \(B\) 的比特 XOR 定义如下，记为 \(A \oplus B\) ：
- 在 \(A \oplus B\) 的二进制表示中，位于位置 \(2^k\) 的数字（ \(k \geq 0\) ）是 \(A \oplus B\) 。(如果 \(A\) 和 \(B\) 的二进制表示中恰好有一个 \(1\) 位于 \(2^k\) 位置（ \(k \geq 0\) ），则该位置的数字为 \(1\) ，否则为 \(0\) 。
例如， \(3 \oplus 5 = 6\) （二进制表示为： \(011 \oplus 101 = 110\) ）。
一般来说， \(k\) 整数 \(p_1, \dots, p_k\) 的 XOR 定义为 \((\cdots ((p_1 \oplus p_2) \oplus p_3) \oplus \cdots \oplus p_k)\) ，可以证明这与 \(p_1, \dots, p_k\) 的阶数无关。\(1 \le N,K \le 2 \times 10^5,1 \le M \le 100\)。
nit 场上给出了正解/bx
根据上一道题的做法，我们发现求的其实是对于每一个 \(p\)
\(\prod\) 表示 \(x\) 的 异或卷积，\(y\) 表示 循环卷积。
然后就和上一道题一样，你对和 FWT 出来系数一定是形如 \(1+y\) 或者 \(1-y\)，容易解出
由于 FWT 是线性变化，所以我们只关心 \([y^0](1+y)^a(1+y)^b\)，再 FWT 回去就可以了。代码。
场上预处理 \((1-y)^b\) 出错也是逆天。
Code
#include <bits/stdc++.h>
using namespace std;
#define ll long long
const int N=2e5+5,M=1<<21,H=998244353;
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;
}
int n,m,k,a[N],lim=1,g[2][N][105],inv2,b[M],ans=0;
struct Nod{
ll x,y;
Nod operator +(const Nod &a) const{return {x+a.x,y+a.y};}
Nod operator -(const Nod &a) const{return {x-a.x,y-a.y};}
Nod operator /(const ll &k) const{return {x/k,y/k};}
}f[M];
void FWT(Nod *a,int lim,int op){
for(int le=1;le<lim;le<<=1)
for(int i=0;i<lim;i+=le<<1)
for(int j=i;j<i+le;j++){
Nod x=a[j],y=a[j+le];
if(op) a[j]=x+y,a[j+le]=x-y;
else a[j]=(x+y)/2,a[j+le]=(x-y)/2;
}
}
void FWT(int *a,int lim,int op){
for(int le=1;le<lim;le<<=1)
for(int i=0;i<lim;i+=le<<1)
for(int j=i;j<i+le;j++){
int x=a[j],y=a[j+le];
if(op) a[j]=adc(x,y),a[j+le]=dec(x,y);
else a[j]=mul(adc(x,y),inv2),a[j+le]=mul(dec(x,y),inv2);
}
}
int main(){
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n>>m>>k;inv2=qpow(2);
for(int i=1;i<=n;i++){
cin>>a[i];
f[0]=f[0]+(Nod){1,0},f[a[i]]=f[a[i]]+(Nod){0,1};
while(lim<=a[i]) lim<<=1;
}
FWT(f,lim,1);
g[0][0][0]=1;
for(int i=1;i<=n;i++)
for(int j=0;j<m;j++)
g[0][i][j]=adc(g[0][i-1][j],g[0][i-1][(j-1+m)%m]);
g[1][0][0]=1;
for(int i=1;i<=n;i++)
for(int j=0;j<m;j++)
g[1][i][j]=dec(g[1][i-1][j],g[1][i-1][(j-1+m)%m]);
for(int i=0;i<lim;i++){
ll nw=f[i].y;
int A=(nw+n)/2,B=(n-nw)/2;
for(int j=0;j<m;j++) add(b[i],mul(g[0][A][j],g[1][B][(m-j)%m]));
}
FWT(b,lim,0);
for(int i=0;i<lim;i++) add(ans,mul(b[i],qpow(i,k)));
cout<<ans;
return 0;
}
浙公网安备 33010602011771号