题意概述
给定$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;
}