P10219 [省选联考 2024] 虫洞

· · 题解

一种不需要高级线性代数知识的低复杂度做法。

我们只考虑 m=0 的情况。实际上，通过之后的分析可以知道 m\neq 0 的情况可以归约为若干个 \sum n 不超过原来的 n 且 m=0 的子问题。

依次加入每种颜色的边。加入第一种颜色时，我们得到了若干个环。

令点 u 所在的环大小为 size(u)。加入第二种颜色时，对于一条边 u\rightarrow v，此时要求 size(u)=size(v)，否则与条件四矛盾。

进一步利用条件四，我们可以刻画第二种颜色的边：将环划分成若干个集合，要求每个集合中所有环大小相等。对于一个集合，令其中的环为 C_1\dots C_p，大小均为 L，令 C_{i,j} 表示环 C_i 中（从任意点开始）顺时针编号的第 j 个点。则 \forall i\in [1,p),j\in [0,L)，存在第二种颜色的边 C_{i,j}\rightarrow C_{i+1,j}。特殊地，\exist d\in [0,L),\forall j\in [0,L)，存在边 C_{p,j}\rightarrow C_{1,(j+d)\bmod L}。

相当于将一个集合中的环合并成了一个“柱”。

再加入第三种颜色。通过一些分析可以得到：将“柱”划分为若干个集合，要求每个集合中所有“柱”同构。连边方式也与之前类似，即为每个“柱”向下一个“柱”的某一种同构方式的对应点连边。

以此类推，我们可以猜测出 k 种颜色的情况：之前的“环”与“柱”都可以被推广为“连通块”，每一轮都是对等价的“连通块”进行合并，形成更大的“连通块”。

实际上，我们可以通过归纳得到一个较为关键的性质：对于一个“连通块”，若看作无标号，那么其中任意两点都是等价的，即存在一种自同构使得这两个点对应，并且这种自同构是唯一的。

因此我们令 f_{i,j} 表示进行 j 步操作能够合并出多少种不同的大小为 i 的无标号“连通块”。有转移：

f_{i,j}\times i\rightarrow f_{i',j+1}

其中要求 i\mid i'。

其意义为这一步合并了 \dfrac{i'}{i} 个两两同构的大小为 i 的“连通块”。因为无标号，所以我们只乘一个 i 表示最后一层到第一层的连边方式。

我们只需要对于每个 i\in [1,n] 计算 f_{i,k}，再进行多项式 exp 即可得到答案。

但 j 这一维可能非常巨大，需要使用一些手段对其进行优化。

将 f 写作生成函数，即令 F_i(x)=\sum f_{i,j}x^j。转移可以改写为：

F_i=\dfrac{1}{1-ix}\sum\limits_{j\mid i}F_j\times jx

进一步将 F_i 写作有理分式，容易发现可以令 F_i 分母为 \prod\limits_{j\mid i}{(1-jx)}。

只需维护 F_i 的分子（可以证明次数不超过分母），最后用线性递推求出 k 次项系数，即可将时间复杂度降为 O(\text{poly}(n,\log k))，使用较好的实现方式可通过本题。

更进一步地，注意到分母形式的特殊性，我们可以将 F_i 写作分式分解形式，即 F_i=\sum\limits_{j\mid i}\dfrac{w_{i,j}}{1-jx}。

按照上述转移式直接暴力维护 w_{i,j}，即可做到 O\left(\sum\dfrac{n}{i}d(i)\right)=O(n\log^2 n)。

这个转移还可以通过类似于“在线高维前缀和”的方式进一步优化，令 s_{i,j} 表示在 i 处只对前 j 个质因子做高维前缀和得到的结果。可以将时间复杂度降为 O(n\log n\log\log n)。

参考代码（O(n\log^2 n)）：

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define pb push_back
const int N=2e6+5,MOD=998244353;
int n,pw[N],id[N],z[N],z1[N];ll m;vector<int> d[N],f[N];
int l,lim,lim1,invL,r[N],inv[N],tmp1[N],tmp2[N],tmp3[N],g[2][N];
void W(int &x,int y) {x+=y;if(x>=MOD) x-=MOD;}
void W(int &x,ll y) {x=(x+y)%MOD;}
int add(int x,int y) {x+=y;return x<MOD?x:x-MOD;}
int qPow(int x,int y)
{int res=1;for(;y;y>>=1,x=1ll*x*x%MOD) if(y&1) res=1ll*res*x%MOD;return res;}
void init(int n)
{
    l=0;lim=1;while(lim<n) ++l,lim*=2;invL=qPow(lim,MOD-2);
    for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<l-1);
    if(lim>lim1)
    {
        for(int i=1;i<lim;++i) inv[i]=i>1?1ll*inv[MOD%i]*(MOD-MOD/i)%MOD:1;
        for(int i=1,t1,t2,t3,t4;i<lim;i*=2)
        {
            t1=qPow(3,(MOD-1)/(i*2));t2=qPow(t1,MOD-2);t3=t4=1;
            for(int j=0;j<i;++j,t3=1ll*t3*t1%MOD,t4=1ll*t4*t2%MOD)
                g[0][i+j]=t3,g[1][i+j]=t4;
        }lim1=lim;
    }
}
void deriv(int n,int a[]) {for(int i=1;i<n;++i) a[i-1]=1ll*a[i]*i%MOD;a[n-1]=0;}
void integ(int n,int a[]) {for(int i=n-1;i;--i) a[i]=1ll*a[i-1]*inv[i]%MOD;a[0]=0;}
void NTT(bool fl,int a[])
{
    for(int i=0;i<lim;++i) if(i<r[i]) swap(a[i],a[r[i]]);
    for(int i=1,t1,t2;i<lim;i*=2) for(int j=0;j<lim;j+=i*2) for(int k=0;k<i;++k)
    {
        t1=a[j+k];t2=1ll*g[fl][i+k]*a[i+j+k]%MOD;
        a[j+k]=add(t1,t2);a[i+j+k]=add(t1,MOD-t2);
    }if(fl) for(int i=0;i<lim;++i) a[i]=1ll*a[i]*invL%MOD; 
}
void polyInv(int n,int a[],int res[])
{
    if(n==1) {res[0]=qPow(a[0],MOD-2);return;}polyInv((n+1)/2,a,res);
    for(int i=0;i<n;++i) tmp1[i]=a[i];for(int i=n;i<lim;++i) tmp1[i]=0;
    init(n*2-1);NTT(0,tmp1);NTT(0,res);
    for(int i=0;i<lim;++i) res[i]=(2-1ll*tmp1[i]*res[i]%MOD+MOD)*res[i]%MOD;
    NTT(1,res);for(int i=n;i<lim;++i) res[i]=0;
}
void polyLn(int n,int a[])
{
    init(n*2-1);for(int i=0;i<lim;++i) tmp2[i]=0;
    polyInv(n,a,tmp2);deriv(n,a);NTT(0,a);NTT(0,tmp2);
    for(int i=0;i<lim;++i) a[i]=1ll*a[i]*tmp2[i]%MOD;
    NTT(1,a);integ(n,a);for(int i=n;i<lim;++i) a[i]=0;
}
void polyExp(int n,int a[],int res[])
{
    if(n==1) {res[0]=1;return;}polyExp((n+1)/2,a,res);
    for(int i=0;i<n;++i) tmp3[i]=res[i];for(int i=n;i<lim;++i) tmp3[i]=0;
    polyLn(n,tmp3);for(int i=0;i<n;++i) tmp3[i]=add(a[i],MOD-tmp3[i]);++tmp3[0];
    NTT(0,tmp3);NTT(0,res);for(int i=0;i<lim;++i) res[i]=1ll*res[i]*tmp3[i]%MOD;
    NTT(1,res);for(int i=n;i<lim;++i) res[i]=0;
}
int main()
{
    scanf("%*d %d %*d %lld",&n,&m);m%=MOD-1;
    for(int i=1;i<=n;++i)
    {
        pw[i]=qPow(i,m);inv[i]=i>1?1ll*inv[MOD%i]*(MOD-MOD/i)%MOD:1;
        for(int j=i;j<=n;j+=i) d[j].pb(i);
    }
    for(int i=1,t;i<=n;++i)
    {
        f[i].resize(d[i].size());if(i==1) {f[i][0]=z[1]=1;continue;}
        for(int j=0;j<d[i].size();++j) id[d[i][j]]=j;
        for(auto j:d[i]) if(j<i) for(int k=0;k<d[j].size();++k)
            W(f[i][id[d[j][k]]],f[j][k]);
        for(int j=0;j<d[i].size();++j)
        {
            if(d[i][j]<i)
            {
                t=1ll*f[i][j]*(MOD-inv[i-d[i][j]])%MOD;
                f[i][j]=t;W(f[i].back(),MOD-t);
            }W(z[i],1ll*f[i][j]*pw[d[i][j]]);
        }z[i]=1ll*z[i]*inv[i]%MOD;for(auto &x:f[i]) x=1ll*x*i%MOD;
    }polyExp(n+1,z,z1);for(int i=1;i<=n;++i) z1[n]=1ll*z1[n]*i%MOD;
    printf("%d\n",z1[n]);return 0;
}