题解：luogu P4948 数列求和

题解：luogu P4948 数列求和

要求：

\[\sum_{i = 1}^{n}{i^k a^i} \]

其中 \(n \leq 10^{18},k \leq 2000\)

这种 \(k\) 次方但是 \(k\) 特别小的一般都是将 \(i^k\) 通过斯特林数展开。

由：

\[x^n=\sum_{i = 0}^{n}{i! \binom{x}{i} {n \brace i}} \]

得：

\[\sum_{i = 1}^{n}{i^k a^i} = \sum_{i = 1}^{n}{a^i \sum_{j = 0}^{k}{j! \binom{i}{j} {k \brace j}}} = \sum_{j = 0}^{k}{j! {k \brace j} \sum_{i = 1}^{n}{a^i {\binom{i}{j}}}} \]

万恶的出题人为了证明这不是多项式题将模数设为了 \(10^9+7\)。

不过斯特林数可以直接 \(O(k^2)\) 求。

考虑后面的怎么求，设 \(s_j = \sum_{i = 1}^{n}{a^i {\binom{i}{j}}}\)，可以得到：

\[s_j = \sum_{i = 1}^{n}{a^i {\binom{i}{j}}} = \sum_{i = 1}^{n}{a^i {\binom{i - 1}{j}}} + \sum_{i = 1}^{n}{a^i \binom{i - 1}{j - 1}} \]

\[\sum_{i = 1}^{n}{a^i {\binom{i - 1}{j}}} = a \sum_{i = 0}^{n - 1}{a^i \binom{i}{j}} = a (s_j - a^n \binom{n}{j} + [j = 0]) \]

\[\sum_{i = 1}^{n}{a^i \binom{i - 1}{j - 1}} = a \sum_{i = 0}^{n - 1}{a^i \binom{i}{j - 1}} = a (s_{j - 1} - a^n \binom{n}{j - 1} + [j = 1]) \]

\[s_j = a (s_j - a^n \binom{n}{j} + [j = 0] + s_{j - 1} - a^n \binom{n}{j - 1} + [j = 1]) \]

\[s_j = \frac{a (-a^n \binom{n}{j} + [j = 0] + s_{j - 1} - a^n \binom{n}{j - 1} + [j = 1])}{1 - a} \]

特别的（其实并不特别）：

\[s_0=\frac{1 - a^{n + 1}}{1 - a} \]

时间复杂度 \(O(k)\)。

观察上式，发现当 \(a = 1\) 时爆掉了，直接除了 \(0\)，那怎么办？

其实直接带入：

\[s_j = a (s_j - a^n \binom{n}{j} + [j = 0] + s_{j - 1} - a^n \binom{n}{j - 1} + [j = 1]) \]

就好了，得到：

\[s_{j - 1} = \binom{n}{j} - [j = 0] + \binom{n}{j - 1} - [j = 1] \\ s_{j} = \binom{n}{j} + \binom{n}{j + 1} - [j = 0] \]

最终答案为：

\[\sum_{i = 0}^{k}{i! {k \brace i} s_i} \]

复杂度 \(O(k^2)\) 或 \(O(k\log k)\)。

一处细节：\(\binom{n}{i} = \binom{n}{i - 1} \times \frac{n - i + 1}{i}\)，这个东西直接递推求就行。

代码：

#include <iostream>using namespace std;#define QED return 0const int mod = 1e9 + 7, N = 2000 + 10;#define int long longint s[N], a, n, k, w[N][N];int qpow(int x, int b)
{int res = 1;while (b){if (b & 1ll) res = res * x % mod;x = x * x % mod;b >>= 1ll;}return res;
}signed main()
{ios :: sync_with_stdio(false), cin.tie(0), cout.tie(0);cin >> n >> a >> k;w[0][0] = 1;for (int i = 1; i <= k; i++){for (int j = 1; j <= k; j++){w[i][j] = (j * w[i - 1][j] % mod + w[i - 1][j - 1]) % mod;}}if (a == 1){int C = 1;for (int i = 0; i <= k; i++){int tmpC = C;if (i != 0) C = C * (n % mod - i + mod + 1) % mod * qpow(i, mod - 2) % mod;tmpC = C;C = C * (n % mod - (i + 1) + mod + 1) % mod * qpow((i + 1), mod - 2) % mod;s[i] = (C + tmpC - (i == 0)) % mod;C = tmpC;}}else{s[0] = (qpow(a, n + 1ll) - a) * qpow(a - 1, mod - 2) % mod;int C = 1;for (int i = 1; i <= k; i++){s[i] = a * (s[i - 1] - C * qpow(a, n) % mod + mod + (i == 1)) % mod;C = C * (n % mod - i + mod + 1) % mod * qpow(i, mod - 2) % mod;s[i] += a * (-C * qpow(a, n) % mod + mod) % mod;s[i] = s[i] * qpow(1 + mod - a, mod - 2) % mod;}}int fac = 1, ans = 0;for (int i = 0; i <= k; i++) fac = fac * max(1ll, i) % mod, ans = (ans + fac * w[k][i] % mod * s[i] % mod) % mod;cout << ans << '\n';QED;
}