给定 \(n, k, p, a_1, a_2, \dots, a_n\)，求 \[ \sum_{1 \leq l \leq r \leq n}\sum_{l \leq i \leq r, 0 \leq t_i \leq a_i} \left(\left(\sum_{i=1}^np^{t_i}t_i^k\right)\left(\sum_{i=1}^n[t_i>0]\right)\right) \]
答案对 \(998244353\) 取模。
\(1 \leq n, k\leq 10^5, 2 \leq p \leq 998244351, 1 \leq a_i \leq 998244351\)。
考虑对固定的 \([l, r]\) 计算答案，它是两个数乘积之和的形式。可以通过维护 \((\sum xy, \sum x, \sum y, \sum 1)\) 四元组来计算答案。
这个玩意是线性变换，写成矩阵形式可以发现，答案最后就是要你求所有区间的矩阵乘积之和。
这个可以直接维护每个位置作为右端点的答案，只要维护当前后缀的四元组即可。
求单个的四元组，可以发现就是要求 \(\sum_{i=0}^{a_j} p^ii^k\)。
对于这个问题，板子 OJ 有一个 单次询问的版本，本题解没有参考上面的做法。
本文主要参考 EI 的评论 以及下面 ecnerwala 的解读。
考虑
\[ \begin{aligned} e^{nx} &= \sum_{i \geq 0}\frac{n^ix^i}{i!}\\ n^k &= [x^k]e^{nx}k!\\ \sum_{i=0}^{n-1}p^ii^k &= [x^k]\sum_{i=0}^{n-1}p^ie^{ix}k!\\ &= [x^k]k!\frac{p^ne^{nx}-1}{pe^x-1}\\ &= [x^k]k!(p^n\frac{e^{nx}-1}{pe^x-1}-\frac{e^{0x}-1}{pe^x-1}) \end{aligned} \]
设 \(S(n)=[x^k]k!\dfrac{e^{nx}-1}{pe^x-1}\)，有 \[ \begin{aligned} p^nS(n) - p^0S(0) &= \sum_{i=0}^{n-1}p^ii^k\\ p^{n+1}S(n+1)-p^nS(n) &= p^nn^k\\ pS(x+1)-S(x) &= x^k \end{aligned} \]
考虑用 \(S(x)\) 来表示 \(S(x+1)\)，即将 \(S(x + 1)\) 在 \(x\) 处展开，得 \[ S(x+1) = \sum_{k \geq 0} \frac{1^k}{k!}{\mathrm D}^kS(x) = e^{\mathrm D}S(x) \] 其中 \(\mathrm D\) 表示微分算子。 则 \[ S(x) = \frac{x^k}{pe^{\mathrm D}-1} \] 这也证明了 \(S(x)\) 是一个 \(k\) 次多项式（考虑生成函数 \(F({\mathrm D}) = \frac{1}{pe^{\mathrm D}-1}\)，高于 \(k\) 次微分 \(x^k\) 就变成 \(0\) 了）。
用求逆算出 \(S(x)\) 后，多点求值即可。
\(O(k\log k + n\log n \log k)\)
void solve() {
int n, k, p;
std::cin >> n >> k >> p;
Poly f(k + 1), g(k + 1);
for (int i = 0; i <= k; i++) f[i] = 1LL * ifac[i] * p % P;
f[0]--;
f = f.inv(k + 1);
for (int i = 0; i <= k; i++) {
g[k - i] = 1LL * fac[k] * ifac[k - i] % P * f[i] % P;
}
Poly pt(n);
for (int &x : pt) std::cin >> x, x++;
f = SegTree(pt).eval(g);
for (int i = 0; i < n; i++) {
f[i] = (1LL * fpow(p, pt[i]) * f[i] + P - g[0]) % P;
}
int ans = 0;
std::array<int, 4> dp = {};
for (int i = 0; i < n; i++) {
inc(dp[0], 1);
auto ndp = dp;
for (auto &x : ndp) x = 1LL * x * pt[i] % P;
ndp[1] = (ndp[1] + 1LL * dp[0] * (pt[i] - 1)) % P;
ndp[2] = (ndp[2] + 1LL * dp[0] * f[i]) % P;
ndp[3] = (ndp[3] + 1LL * (dp[0] + dp[1]) * f[i] + 1LL * dp[2] * (pt[i] - 1)) % P;
dp = ndp;
inc(ans, dp[3]);
}
std::cout << ans << "\n";
}最后修改于 2022-05-12