【学习笔记】[ARC145F] Modulo Sum of Increasing Sequences

经典题，简单记录一下。

单位根反演好题。

提示：是照搬的 这篇题解 的做法，只是加了一点小小的解释。

首先，做等价变换：给第$i$个位置加上$i-1$，问题转化为了求单调递增序列，即从$[0,N+M-1]$中选$N$个不同的数，使得这些数模$\operatorname{mod}$的值为$k$。

这事实上是一个二元$GF$的问题。先不考虑选$N$个数的限制，记$n=N+M$，写出答案的生成函数形式：$\prod_{i=0}^{n-1}(1+x^i)$

记$k=\operatorname{mod}$，那么我们要求所有$ik+r$项的系数和，这通常可以用单位根反演来解决。

对于任意多项式$f(x)$，我们有：

$\begin{aligned} \sum_{k|i}[x^i]f(x)&=\sum_{i=0}^{n-1}[x^i]f(x)\frac{1}{k}\sum_{j=0}^{k-1}w_k^{ij}\\&=\frac{1}{k}\sum_{i=0}^{n-1}a_i\sum_{j=0}^{k-1}(w_k^j)^i\\&=\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^{n-1}a_i(w_k^j)^i\\&=\frac{1}{k}\sum_{j=0}^{k-1}f(w_k^j) \end{aligned}$

类似的：$\sum[x^{ik+r}]f(x)=\frac{1}{k}\sum_{j=0}^{k-1}w_k^{-jr}f(w_k^j)$

注意这里的$k$不是$2$的次幂，因此不能直接计算。要想办法把单位根搞掉。

注意到$x^n-1=\prod_{i=0}^{n-1}(x-w_n^i)$，令$x=-1$，那么$\prod_{i=0}^{n-1}(w_n^i+1)=1-(-1)^n$。

考虑$k|n$的情况。 记$d=\gcd(j,k)$，此时$\gcd(\frac{j}{d},\frac{k}{d})=1$：

$$ \begin{aligned} f(w_k^j)&=\prod_{i=0}^{n-1}(1+w_k^{ij})\\&=\prod_{i=0}^{\frac{k}{d}-1}(1+w_{\frac{k}{d}}^{i})^{\frac{nd}{k}}\\&=(1-(-1)^{\frac{k}{d}})^{\frac{nd}{k}} \end{aligned}$$

带入原式即：$\frac{1}{k}\sum_{j=0}^{k-1}w_k^{-jr}(1-(-1)^{\frac{k}{d}})^{\frac{nd}{k}}$，现在考虑怎么把前面的单位根搞掉。看到$\gcd$考虑枚举$d$，即：

$$LHS=\frac{1}{k}\sum_{d|k}(1-(-1)^{\frac{k}{d}})^{\frac{nd}{k}}\sum_{\gcd(j,k)=d}w_k^{-jr}$$

后面那一坨看起来就很好化简：

$$ \begin{aligned} \sum_{\gcd(j,k)=d}w_k^{-jr}&=\sum_{i=1}^{\frac{k}{d}}w_k^{-idr}[\gcd(i,\frac{k}{d})=1]\\&=\sum_{i=1}^{\frac{k}{d}}w_k^{-idr}\sum_{j|\gcd(i,\frac{k}{d})}\mu(j)\\&=\sum_{j|\frac{k}{d}}\mu(j)\sum_{i=1}^{\frac{k}{dj}}w_k^{-ijdr}\\&=\sum_{j|\frac{k}{d}}\mu(j)\sum_{i=1}^{\frac{k}{dj}}w_\frac{k}{dj}^{-ri}\\&=\sum_{j|\frac{k}{d}}\mu(j)[\frac{k}{dj}|r] \end{aligned} $$

最后一步是 逆用单位根反演 。因此最终答案为：

$$\begin{aligned} LHS=\frac{1}{k}\sum_{d|k}(1-(-1)^{\frac{k}{d}})^{\frac{nd}{k}}g(d,k) \end{aligned}$$

其中$g(d,k)$是固定的系数。现在考虑二元$GF$，即$[x^{ik+r}y^n]\prod_{i=0}^{n-1}(1+x^iy)$。前面的推导都一模一样（把$y$当成参量来处理），不难验证答案为：

$$LHS=\frac{1}{k}\sum_{d|k}(1-(-y)^{\frac{k}{d}})^{\frac{nd}{k}}g(d,k)$$

对于$n\nmid k$的情况，可以将剩余的不足$k$项暴力背包，即设$f_{i,j}$表示选了$i$个数，并且模$k$等于$j$的方案数。可以在$O(k^3)$的时间内处理出来，最后合并答案也是$O(k^3)$的。

#include<bits/stdc++.h>
#define ll long long
#define pb push_back
#define fi first
#define se second
using namespace std;
const int mod=998244353;
const int N=505;
const int M=2e6+5;
int n,m,p,n2,mx;
ll now[N][N];
ll fac[M],inv[M],f[N],g[N][N],res[N];
ll fpow(ll x,ll y=mod-2){
    ll z(1);
    for(;y;y>>=1){
        if(y&1)z=z*x%mod;
        x=x*x%mod;
    }return z;
}
void init(int n){
    fac[0]=1;for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;
    inv[n]=fpow(fac[n]);for(int i=n;i>=1;i--)inv[i-1]=inv[i]*i%mod;
}
int mu[N],vs[N];
int pr[N],cnt;
void init2(int n){
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(vs[i]==0)pr[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&pr[j]<=n/i;j++){
            vs[i*pr[j]]=1;
            if(i%pr[j]==0){
                mu[i*pr[j]]=0;
                break;
            }mu[i*pr[j]]=-mu[i];
        }
    }
}
ll binom(int x,int y){
    if(x<0||y<0||x<y)return 0;
    return fac[x]*inv[y]%mod*inv[x-y]%mod;
}
void add(ll &x,ll y){
    x=(x+y)%mod;
}
vector<int>v;
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    cin>>n>>m>>p,n2=n+m;
    init(n2),init2(p);
    mx=p*(n2/p);
    for(int i=1;i<=p;i++){
        if(p%i==0)v.pb(i);
    }
    for(int k=0;k<p;k++){
        for(auto d:v){
            for(auto j:v){
                if(p/d%j==0&&k%(p/d/j)==0){
                    add(g[k][d],mu[j]*p/d/j);
                }
            }
        }
    }
    now[0][0]=1;
    for(int i=mx;i<n2;i++){
        int s=i%p;
        for(int j=min(n2-mx,n);j>=1;j--){
            for(int k=0;k<p;k++){
                add(now[j][k],now[j-1][(k-s+p)%p]);
            }
        }
    }
    for(int i=0;i<=min(n,n2-mx);i++){
        memset(f,0,sizeof f);
        for(auto d:v){
            if((n-i)%(p/d))continue;
            int t=(n-i)/(p/d);
            for(int k=0;k<p;k++){
                if(p/d%2==1||t%2==0){
                    add(f[k],binom(d*mx/p,t)*g[k][d]);
                }
                else{
                    add(f[k],-binom(d*mx/p,t)*g[k][d]);
                }
            }
        }
        for(int j=0;j<p;j++){
            for(int k=0;k<p;k++){
                add(res[(j+k)%p],f[j]*now[i][k]);
            }
        }
    }for(int i=0;i<p;i++){
        res[i]=res[i]*fpow(p)%mod;
    }int tmp=(ll)n*(n-1)/2%p;
    for(int i=0;i<p;i++){
        cout<<(res[(i+tmp)%p]+mod)%mod<<" ";
    }
}