序列求和 V4

题意概述

给定$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;
}