Codeforces 622F 拉格朗日插值 + 线性筛

题意

传送门 Codeforces 622F The Sum of the k-th Powers

题解

$\sum _{i=1}^n{i}^k$ 是一个关于 $n$ 的 $k+1$ 次多项式。用 $1,2,3\cdots k+2$ 共 $k+2$ 个值，进行拉格朗日插值。令 $f(n)=\sum _{i=1}^n{i}^k$。得到
$$f(x)=\sum _{i=1}^{k+2}f(i)\prod _{j=1,j\ne i}^{k+2}\frac{x-j}{i-j}$$ 分子可以表示为
$$\frac{\prod _{j=1}^{k+2}x-j}{x-i}$$ 可以通过预处理前缀与后缀积 $O(1)$ 得到。分母是
$$(i-1)!(k+2-i)!(-1{)}^{k+2-i}$$ 可以通过预处理阶乘的逆元 $O(1)$ 得到。对于 $f(i)$，满足 $f(i)=f(p)f(i/p),p|i$，可以通过线性筛递推。由于 $[1,n]$ 内质数规模为 $O(n/\ln n)$，且仅当筛到质数时进行快速幂，故求解 $f(i)$ 总时间复杂度接近于线性。

总时间复杂度 $O(k)$。

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
constexpr int MOD = (int)1e9 + 7;

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n, k;
    cin >> n >> k;
    int mx = k + 2;
    vector<ll> fac(mx + 1), inv(mx + 1), invf(mx + 1);
    fac[0] = fac[1] = 1;
    inv[1] = 1;
    invf[0] = invf[1] = 1;
    for (int i = 2; i <= mx; ++i) {
        fac[i] = fac[i - 1] * i % MOD;
        inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
        invf[i] = invf[i - 1] * inv[i] % MOD;
    }

    vector<ll> pre(mx + 2), suf(mx + 2);
    pre[0] = suf[mx + 1] = 1;
    for (int i = 1; i <= mx; ++i) {
        pre[i] = pre[i - 1] * (n - i) % MOD;
    }
    for (int i = mx; i > 0; --i) {
        suf[i] = suf[i + 1] * (n - i) % MOD;
    }

    auto power = [&](ll x, int n) {
        ll res = 1;
        while (n > 0) {
            if (n & 1) (res *= x) %= MOD;
            (x *= x) %= MOD, n >>= 1;
        }
        return res;
    };

    vector<ll> f(mx + 1);
    vector<int> prime, minp(mx + 1);
    f[1] = 1;
    for (int i = 2; i <= mx; ++i) {
        if (minp[i] == 0) {
            minp[i] = i;
            f[i] = power(i, k);
        }
        for (int j = 0, x; j < (int)prime.size() && (x = prime[j] * i) <= mx; ++j) {
            minp[x] = prime[j];
            f[x] = f[i] * f[prime[j]] % MOD;
            if (minp[i] == prime[j]) break;
        }
    }
    for (int i = 2; i <= mx; ++i) (f[i] += f[i - 1]) %= MOD;

    ll res = 0;
    for (int i = 1; i <= mx; ++i) {
        ll up = pre[i - 1] * suf[i + 1] % MOD;
        ll down = invf[i - 1] * invf[k + 2 - i] % MOD;
        if ((k + 2 - i) & 1) down = MOD - down;
        (res += up * down % MOD * f[i] % MOD) %= MOD;
    }

    res = (res + MOD) % MOD;
    cout << res << '\n';

    return 0;
}