ARC145F Modulo Sum of Increasing Sequences
Ender32k
2023-07-19 21:18:07
[更好的阅读体验。](https://www.cnblogs.com/Ender32k/p/17569358.html)
为数不多不用多项式科技的单位根反演题。
$A$ 不降比较难搞，所以首先令 $B_i=A_i+i-1$，则 $B$ 单调递增。转化为对任意的 $k\in [0,\text{MOD}-1]$，求在 $[0,N+M-1]$ 中选 $N$ 个**不同**的数，总和对 $\text{MOD}$ 取模为 $k$ 的方案数。记 $p=\text{MOD},n=N+M$。
列出选数的 OGF：
$$F(x)=\prod\limits_{i=0}^{n-1}(1+x^i)$$
我们现在不考虑长度为 $N$ 的限制，那么答案就是：
$$\text{ans}=\sum\limits_{p\mid i-k}[x^i]F(x)$$
根据单位根反演，对于一个 $n$ 次多项式 $f(x)$，其所有 $k$ 的倍数次方的系数之和：
$$\begin{aligned}\sum\limits_{i=0}^{\lfloor\frac{n}{k}\rfloor}[x^{ik}]f(x)&=\sum\limits_{i=0}^n[k\mid i][x^i]f(x)\\&=\sum\limits_{i=0}^n[x^i]f(x)\frac{1}{k}\sum\limits_{j=0}^{k-1}\omega_{k}^{ij}\\&=\frac{1}{k}\sum\limits_{i=0}^na_i\sum\limits_{j=0}^{k-1}\omega_{k}^{ij}\\&=\frac{1}{k}\sum\limits_{i=0}^n\sum\limits_{j=0}^{k-1}a_i(\omega_k^j)^i\\&=\frac{1}{k}\sum\limits_{j=0}^{k-1}f(\omega_k^j)\end{aligned}$$
同理，当我们要求 $ik+r$ 的系数：
$$\sum\limits_{i=0}^{\lfloor\frac{n}{k}\rfloor}[x^{ik+r}]f(x)=\frac{1}{k}\sum\limits_{j=0}^{k-1}f(\omega_k^j)\omega_{k}^{-jr}$$
于是对 $\text{ans}$ 使用单位根反演：
$$\begin{aligned}\text{ans}&=\sum\limits_{p\mid i-k}[x^i]F(x)\\&=\frac{1}{p}\sum\limits_{i=0}^{p-1}F(\omega_p^i)\omega_p^{-ik}\\&=\frac{1}{p}\sum\limits_{i=0}^{p-1}\omega_p^{-ik}\prod\limits_{j=0}^{n-1}(1+\omega_p^{ij})\end{aligned}$$
考虑 $p\mid n$ 的情况，记 $d=\gcd(i,p)$，那么 $\omega_p^i$ 在 $0$ 到 $n-1$ 次方有纯循环节 $p$，而 $\omega_{p}^i$ 有循环节 $d$，并考虑 $\omega_p^{ij}=\omega_{\frac{p}{d}}^{\frac{ij}{d}}$，于是枚举 $\frac{i}{d}$：
$$\begin{aligned}\prod\limits_{j=0}^{n-1}(1+\omega_p^{ij})&=\left(\prod\limits_{j=0}^{p-1}(1+\omega_p^{ij})\right)^{\frac{n}{p}}\\&=\left(\prod\limits_{j=0}^{\frac{p}{d}-1}(1+\omega_\frac{p}{d}^{ij})^d\right)^{\frac{n}{p}}\\&=\left(\prod\limits_{j=0}^{\frac{p}{d}-1}(1+\omega_{\frac{p}{d}}^{ij})\right)^{\frac{dn}{p}}\end{aligned}$$
然后考虑**分圆多项式**：
$$x^n-1=\prod\limits_{i=0}^{n-1}(x-\omega_n^i)$$
代入 $x=-1$ 得：
$$\begin{aligned}(-1)^n-1&=\prod\limits_{i=0}^{n-1}(-1-\omega_n^i)\\(-1)^n-1&=(-1)^{n}\prod\limits_{i=0}^{n-1}(1+\omega_n^i)\\1-(-1)^n&=\prod\limits_{i=0}^{n-1}(1+\omega_n^i)\end{aligned}$$
然后将右式代入上面的式子：
$$\begin{aligned}\prod\limits_{j=0}^{n-1}(1+\omega_p^{ij})&=\left(\prod\limits_{j=0}^{\frac{p}{d}-1}(1+\omega_{\frac{p}{d}}^{ij})\right)^{\frac{dn}{p}}\\&=(1-(-1)^\frac{p}{d})^{\frac{dn}{p}}\end{aligned}$$
然后代入 $\text{ans}$：
$$\begin{aligned}\text{ans}&=\frac{1}{p}\sum\limits_{i=0}^{p-1}\omega_p^{-ik}\prod\limits_{j=0}^{n-1}(1+\omega_p^{ij})\\&=\frac{1}{p}\sum\limits_{i=0}^{p-1}\omega_p^{-ik}(1-(-1)^{\frac{p}{d}})^{\frac{dn}{p}}\end{aligned}$$
这样我们就消掉了一些单位根，接着我们考虑 $\omega_p^{-ik}$ 如何处理。由于 $d=\gcd(p,i)$，我们考虑枚举 $d$：
$$\begin{aligned}\text{ans}&=\frac{1}{p}\sum\limits_{i=0}^{p-1}\omega_p^{-ik}(1-(-1)^{\frac{p}{d}})^{\frac{dn}{p}}\\&=\frac{1}{p}\sum\limits_{i=0}^{p-1}(1-(-1)^{\frac{p}{d}})^{\frac{dn}{p}}\sum\limits_{\gcd(i,p)=d}\omega_p^{-ik}\end{aligned}$$
把后面那个求和单独拎出来设为 $g(d)$，考虑莫反以及 $\omega_n^k=\omega_{\alpha n}^{\alpha k}$：
$$\begin{aligned}\sum\limits_{\gcd(i,p)=d}\omega_p^{-ik}&=\sum\limits_{i=0}^{p-1}[\gcd(i,p)=d]\omega_p^{-ik}\\&=\sum\limits_{i=1}^{\frac{p}{d}-1}[\gcd(i,\frac{p}{d})=1]\omega_{\frac{p}{d}}^{-ik}\\&=\sum\limits_{i=1}^{\frac{p}{d}-1}\sum\limits_{j\mid \frac{p}{d},j\mid i}\mu(j)\omega_{\frac{p}{d}}^{-ik}\\&=\sum\limits_{j\mid \frac{p}{d}}\mu(j)\sum\limits_{i=0}^{\frac{p}{dj}-1}\omega_{\frac{p}{dj}}^{-ik}\\&=\sum\limits_{j\mid \frac{p}{d}}\mu(j)[\frac{p}{dj}\mid k]\frac{p}{dj}\end{aligned}$$
最后那一步是逆用单位根反演。
于是我们可以对每个 $d$，$O(p)$ 预处理出 $g(d)$ 的值，总共是 $O(p^2)$ 的。
最后再代入答案式子：
$$\begin{aligned}\text{ans}&=\frac{1}{p}\sum\limits_{i=0}^{p-1}(1-(-1)^{\frac{p}{d}})^{\frac{dn}{p}}\sum\limits_{\gcd(i,p)=d}\omega_p^{-ik}\\&=\frac{1}{p}\sum\limits_{i=0}^{p-1}(1-(-1)^{\frac{p}{d}})^{\frac{dn}{p}}g(d)\end{aligned}$$
根据**二项式定理**，上面的式子是容易求的，展开即可。
这是不考虑 $N$ 长度限制的情况，现在我们加上长度限制，那么相当于给 OGF 多加一个元：
$$F(x,y)=\prod\limits_{i=0}^{n-1}(1+x^iy)$$
于是根据上面的推导，我们要求的就是：
$$[y^N]\frac{1}{p}\sum\limits_{i=0}^{p-1}(1-(-y)^{\frac{p}{d}})^{\frac{dn}{p}}g(d)$$
于是 $p\mid n$ 时，我们可以 $O(p^2)$ 解决上述问题。
考虑 $p\nmid n$ 的情况。
显然我们可以先 $O(p^2)$ 解决只能取 $0$ 到 $p\lfloor\frac{n}{p}\rfloor$ 的数的情况，接下来剩下 $n-p\lfloor\frac{n}{p}\rfloor<p$ 个数，于是对于剩下的数我们可以 $O(p^3)$ 暴力 dp 出所有的方案数。具体地，令 $f_{i,j,k}$ 表示从 $p\lfloor\frac{n}{p}\rfloor+1$ 考虑到 $i$，取了 $j$ 个数，和对 $p$ 取模为 $k$ 的方案数，暴力转移，最后暴力合并答案即可。可能需要滚动数组之类的东西。
复杂度 $O(N+M+p^3)$。
```cpp
// Problem: F - Modulo Sum of Increasing Sequences
// Contest: AtCoder - AtCoder Regular Contest 145
// URL: https://atcoder.jp/contests/arc145/tasks/arc145_f
// Memory Limit: 1024 MB
// Time Limit: 4000 ms
//
// Powered by CP Editor (https://cpeditor.org)
#include <bits/stdc++.h>
#define int long long
using namespace std;
namespace vbzIO {
char ibuf[(1 << 20) + 1], *iS, *iT;
#if ONLINE_JUDGE
#define gh() (iS == iT ? iT = (iS = ibuf) + fread(ibuf, 1, (1 << 20) + 1, stdin), (iS == iT ? EOF : *iS++) : *iS++)
#else
#define gh() getchar()
#endif
#define mt make_tuple
#define mp make_pair
#define fi first
#define se second
#define pc putchar
#define pb emplace_back
#define ins insert
#define era erase
typedef tuple<int, int, int> tu3;
typedef pair<int, int> pi;
inline int rd() {
char ch = gh();
int x = 0;
bool t = 0;
while (ch < '0' || ch > '9') t |= ch == '-', ch = gh();
while (ch >= '0' && ch <= '9') x = (x << 1) + (x << 3) + (ch ^ 48), ch = gh();
return t ? ~(x - 1) : x;
}
inline void wr(int x) {
if (x < 0) x = ~(x - 1), putchar('-');
if (x > 9) wr(x / 10);
putchar(x % 10 + '0');
}
}
using namespace vbzIO;
const int M = 1e3 + 100;
const int N = 2e6 + 200;
const int P = 998244353;
int n, m, p, mx, fac[N], inv[N], ifac[N], f[M][M], vl[M][M], cf[M];
int ct, mu[M], pr[M], vs[M], ans[M];
vector<int> di[M];
void init(int lim) {
fac[0] = ifac[0] = inv[1] = 1;
for (int i = 1; i <= lim; i++) {
if (i > 1) inv[i] = 1ll * inv[P % i] * (P - P / i) % P;
fac[i] = 1ll * fac[i - 1] * i % P;
ifac[i] = 1ll * ifac[i - 1] * inv[i] % P;
}
}
int C(int n, int m) { return (n >= m && m >= 0) ? (1ll * fac[n] * ifac[m] % P * ifac[n - m] % P) : 0; }
void gtmu(int lim) {
mu[1] = 1;
for (int i = 2; i <= lim; i++) {
if (!vs[i]) mu[pr[++ct] = i] = P - 1;
for (int j = 1; j <= ct && i * pr[j] <= lim; j++) {
vs[i * pr[j]] = 1;
if (i % pr[j] == 0) {
mu[i * pr[j]] = 0;
break;
}
mu[i * pr[j]] = P - mu[i];
}
}
}
signed main() {
n = rd(), m = rd(), p = rd(), init(N - 10), gtmu(M - 10);
for (int i = 1; i <= p; i++)
for (int j = i; j <= p; j += i)
di[j].pb(i);
for (int d : di[p]) {
for (int k = 0; k < p; k++)
for (int j : di[p / d])
if (d * j * k % p == 0)
(vl[d][k] += 1ll * p / d / j * mu[j] % P) %= P;
}
n += m, m = n - m, mx = p * (n / p), f[0][0] = 1;
for (int i = 1; i <= n - mx; i++)
for (int j = i; j; j--)
for (int k = 0; k < p; k++)
(f[j][k] += f[j - 1][(k - (i + mx - 1) % p + p) % p]) %= P;
for (int k = 0; k <= min(n - mx, m); k++) {
int res = m - k;
memset(cf, 0, sizeof(int) * (p + 10));
for (int d : di[p]) {
int b = p / d, c;
if (res % b) continue;
if (b & 1) c = inv[p];
else c = ((res / b) & 1) ? P - inv[p] : inv[p];
c = 1ll * c * C(d * mx / p, res / b) % P;
for (int i = 0; i < p; i++)
(cf[i] += 1ll * c * vl[d][i] % P) %= P;
}
for (int i = 0; i < p; i++)
for (int j = 0; j < p; j++)
(ans[(i + j) % p] += 1ll * cf[i] * f[k][j] % P) %= P;
}
int _ = 1ll * m * (m - 1) / 2 % p;
for (int i = 0; i < p; i++)
wr(ans[(i + _) % p]), pc(' ');
return 0;
}
```