题意概述
给定$n$和$k$,求$\sum_{i=1}^n i^k \bmod 10^9+7$。有$T$组数据。
数据范围:$1 \le T \le 500, \; 1 \le n \le 10^{18}, \; 1 \le k \le 50000$。
算法分析1
设答案为$f(n)$。这是一个$(k+1)$次多项式,只要确定$(k+2)$个点就可以确定$f$。
我们可以令$x=1\ldots k+2$,在$O(k\log k)$时间内计算出它们对应的$y$。
利用拉格朗日插值法。假设有$n$个点$(x_i, y_i)$,那么
$$f(x)=\sum_{i=1}^n {\prod_{j \in [1, n] \land j \neq i} (x-x_j) \over \prod_{j \in [1, n] \land j \neq i} (x_i-x_j)}y_i$$
考虑分子部分。在给定$x$的情况下,这一部分可以$O(k)$预处理$O(1)$求值。
考虑分母部分。它形如$\prod_{j \in [i-k-2, i-1] \land j \neq 0} j$,可以在插值时$O(1)$维护。
总时间复杂度为$O(Tk\log k)$。
代码1
#include <cstdio> #define int long long void read(int &n) { char c; while ((c = getchar()) < '0' || c > '9') ; n = c - '0'; while ((c = getchar()) >= '0' && c <= '9') (n *= 10) += c - '0'; } static const int K = 50005; static const int MOD = 1000000007; int n, k, a[K], p[K], q[K], inv[K], base[K]; int power(int a, int b) { int ret = 1; a %= MOD, b %= MOD - 1; while (b) b & 1 && ((ret *= a) %= MOD), (a *= a) %= MOD, b >>= 1; return ret; } int calc(int n) { if (n <= k + 2) return a[n]; n %= MOD; int w = power(base[k + 2], MOD - 2), ans = 0; p[0] = q[k + 3] = 1; for (int i = 1; i <= k + 2; ++ i) p[i] = p[i - 1] * (n - i) % MOD; for (int i = k + 2; i; -- i) q[i] = q[i + 1] * (n - i) % MOD; for (int i = 1; i <= k + 2; ++ i) (ans += a[i] * w % MOD * p[i - 1] % MOD * q[i + 1]) %= MOD, w = w * (i - k - 2) % MOD * inv[i] % MOD; return (ans + MOD) % MOD; } signed main() { int T; read(T), base[1] = 1; for (int i = 1; i < K; ++ i) inv[i] = power(i, MOD - 2); for (int i = 2; i < K; ++ i) base[i] = base[i - 1] * (MOD + 1 - i) % MOD; while (T --) { read(n), read(k); for (int i = 1; i < k + 3; ++ i) a[i] = (a[i - 1] + power(i, k)) % MOD; printf("%lld\n", calc(n)); } return 0; }
算法分析2
此题也可以用伯努利数解决,因为
$$\sum_{i=1}^n i^k={1 \over k+1} \sum_{i=0}^k (-1)^i {k+1 \choose i} B_in^{k+1-i}$$
其中伯努利数由其指数型生成函数${x \over e^x-1}$定义。易知
$$\sum_{i=0}^{\infty} B_i{x^i \over i!}={x \over e^x-1}={x \over \sum_{i=1}^{\infty} {x^i \over i!}}={1 \over \sum_{i=0}^{\infty} {x^i \over (i+1)!}}$$
只要对分母求多项式逆元,即可得到伯努利数,时间复杂度$O(k\log k)$。其它部分均可$O(k)$预处理,答案也可以$O(k)$计算。因此总时间复杂度是$O(k\log k+Tk)$。
代码2
#include <algorithm> #include <cmath> #include <cstdio> #include <cstring> #define int long long static int const N = 200000; static int const MOD = 1000000007; static int const M = 31622; static double const PI = acos(-1); int b[N], rev[N], fac[N], inv[N]; class complex { private: double x, y; public: complex(double _x = 0, double _y = 0) : x(_x), y(_y) {} double real() { return x; } friend complex operator+(complex const &a, complex const &b) { return complex(a.x + b.x, a.y + b.y); } friend complex operator-(complex const &a, complex const &b) { return complex(a.x - b.x, a.y - b.y); } friend complex operator*(complex const &a, complex const &b) { return complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); } } wn[N], A[N], B[N], C[N], D[N], E[N], F[N], G[N]; int power(int a, int b) { int ret = 1; for (; b; b >>= 1) b & 1 && ((ret *= a) %= MOD), (a *= a) %= MOD; return ret; } int init(int m) { int n = 1, l = 0; for (; n <= m; n <<= 1, ++l) ; for (int i = 1; i < n; ++i) rev[i] = rev[i >> 1] >> 1 | (i & 1) << l - 1; for (int i = 0; i < n >> 1; ++i) wn[i] = complex(cos(2 * PI / n * i), sin(2 * PI / n * i)); return n; } void fft(complex *a, int n, bool inv = 0) { for (int i = 0; i < n; ++i) if (i < rev[i]) std::swap(a[i], a[rev[i]]); for (int i = 1; i < n; i <<= 1) for (int j = 0; j < n; j += i << 1) for (int k = 0; k < i; ++k) { complex x = a[j + k], y = wn[n / (i << 1) * k] * a[j + k + i]; a[j + k] = x + y, a[j + k + i] = x - y; } if (inv) { std::reverse(a + 1, a + n); for (int i = 0; i < n; ++i) a[i] = a[i].real() / n; } } int get_mod(double t) { return (int)(fmod(t, MOD) + MOD + 0.5) % MOD; } void get_inv(int *f, int *g, int n) { if (n == 1) return void(g[0] = power(f[0], MOD - 2)); get_inv(f, g, n + 1 >> 1); int len = init(n << 1); for (int i = 0; i < n + 1 >> 1; ++i) A[i] = g[i] / M, B[i] = g[i] % M; for (int i = n + 1 >> 1; i < len; ++i) A[i] = B[i] = 0; for (int i = 0; i < n; ++i) C[i] = f[i] / M, D[i] = f[i] % M; for (int i = n; i < len; ++i) C[i] = D[i] = 0; fft(A, len), fft(B, len), fft(C, len), fft(D, len); for (int i = 0; i < len; ++i) { E[i] = 0 - A[i] * C[i]; F[i] = 0 - A[i] * D[i] - B[i] * C[i]; G[i] = 2 - B[i] * D[i]; } fft(E, len, 1), fft(F, len, 1), fft(G, len, 1); for (int i = 0; i < n; ++i) { int x = get_mod(E[i].real()) * M % MOD * M % MOD; int y = get_mod(F[i].real()) * M % MOD; int z = get_mod(G[i].real()); int w = (x + y + z) % MOD; C[i] = w / M, D[i] = w % M; } for (int i = n; i < len; ++i) C[i] = D[i] = 0; fft(C, len), fft(D, len); for (int i = 0; i < len; ++i) { E[i] = A[i] * C[i]; F[i] = A[i] * D[i] + B[i] * C[i]; G[i] = B[i] * D[i]; } fft(E, len, 1), fft(F, len, 1), fft(G, len, 1); for (int i = 0; i < n; ++i) { int x = get_mod(E[i].real()) * M % MOD * M % MOD; int y = get_mod(F[i].real()) * M % MOD; int z = get_mod(G[i].real()); g[i] = (x + y + z) % MOD; } } int get_c(int n, int m) { return fac[n] * inv[m] % MOD * inv[n - m] % MOD; } signed main() { int T; fac[0] = 1; for (int i = 1; i < N; ++i) fac[i] = fac[i - 1] * i % MOD; inv[N - 1] = power(fac[N - 1], MOD - 2); for (int i = N - 1; i; --i) inv[i - 1] = inv[i] * i % MOD; get_inv(inv + 1, b, 50001); for (int i = 0; i < 50001; ++i) (b[i] *= fac[i]) %= MOD; for (scanf("%lld", &T); T--;) { int n, k, ans = 0; scanf("%lld%lld", &n, &k), n %= MOD; for (int i = k, f = k & 1 ? -1 : 1, p = n; ~i; --i, f = -f, (p *= n) %= MOD) ans += (MOD + f) * get_c(k + 1, i) % MOD * b[i] % MOD * p % MOD; printf("%lld\n", ans % MOD * power(k + 1, MOD - 2) % MOD); } return 0; }