拉格朗日插值法

拉格朗日插值法是一种多项式插值方法。简单来说，给定 $n$ 个点，它可以唯一地构造出经过这些点的 $n-1$ 次函数。
如何实现呢？其实非常地暴力：记我们给定的点为 $(a_i,b_i)$，我们只需要构造出一些函数 $f_i(x)$ 使得 $f_i(x)$ 满足 $f(a_i)=b_i$，$f(a_j)=0$，$j̸=i$ 即可。
接下来介绍“开关函数” $p_i(x)$。$p_i(x)$ 满足当 $x=a_i$ 时 $p_i(x)=1$，$x=a_j$，$j̸=i$ 时 $p_i(x)=0$。这样 $f_i(x)=p_i(x)b_i$。
如何实现上述功能呢？我们考虑一个简单的数学性质：$i/i=0$，$0/i=0$，$i̸=0$。这样我们只需要让 $x=a_i$ 时分子分母一定相等，$x$ 为其他的 $a$ 值时分母一定出现 $0$ 即可。可以得到：$p_i={\prod }_{j=1,j̸=i}^n\frac{x-a_j}{a_i-a_j}$。考虑这样的 $p_i$ 为什么拥有上述性质：当 $x=a_i$ 时，无论 $j$ 为何值，分子分母一定相等，即 $p_i=1$；当 $x=a_m$，$m̸=i$ 时，由于 $j̸=i$，故必然存在 $j=m$，此时分子为 $0$，$p_i=0$。
综上所述，我们可以给出一般地过点集 { ( a i , b i ) } \{(a_i,b_i)\} {(ai​,bi​)} 的多项式插值公式：
$$\begin{array}{lll}F(n)& =& {\sum }_{i=1}^n{f}_i(n)\\ & =& {\sum }_{i=1}^n{b}_i{\prod }_{j=1,j̸=i}^n\frac{n-a_j}{a_i-a_j}\end{array}$$
注意其中 $i$，$j$ 和 $n$ 的区别。

例题

The sum of the K-th Powers

给出两个非负整数 $n$，$k$，求
$$\sum _{i=1}^n{i}^k$$
答案对 $10^9+7$ 取模。
$n\le 10^9$，$k\le 10^6$。

题解

暴力计算显然实现复杂度无法承受。
记 $F(n)={\sum }_{i=1}^n{i}^k$。
首先分析当 $k$ 较小的情况。$k=1$ 时，$F(n)=n(n+1)/2$；$k=2$ 时，$F(n)=n(n+1)(2n+1)/6$；依此类推，可以发现 $F(n)$ 是关于 $n$ 的 $k+1$ 次多项式。
既然如此，我们可以通过先暴力算出一部分多项式上的点，然后通过插值得出多项式关于 $n$ 的函数形式，再直接通过函数求值。
我们知道，确定一个 $n$ 次函数需要 $n+1$ 个点，故我们首先暴力算出 $F(1),F(2),...,F(k+2)$。然后，我们将这些点代入拉格朗日插值公式得
$$F(n)=\sum _{i=1}^{k+2}F(i)\prod _{j=1,j̸=i}^{k+2}\frac{n-j}{i-j}$$
然后就到了愉快的拆式子时间了。
$$F(n)=\sum _{i=1}^{k+2}F(i)\frac{({\prod }_{j=1}^{k+2}(n-j))/(n-i)}{(i-1)!(k+2-i)!(-1{)}^{k+2-i}}$$
这样只需要通过预处理得出分母中阶乘的逆元，并通过费马小定理求出 $i-1$ 的逆元即可在 $O(k\log k)$ 的时间内得出答案。

代码

#include <iostream>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const unsigned MOD = 1E9 + 7;
const unsigned MAXN = 1E6;
const unsigned N = MAXN + 10;
unsigned n, k;
unsigned inv[N];
unsigned f[N];
void init() {
    inv[0] = inv[1] = 1;
    for (int i = 2; i <= k + 2; i++) {
        inv[i] = (-ll(MOD / i) * inv[MOD % i] % MOD + MOD) % MOD;
    }
    for (int i = 2; i <= k + 2; i++) {
        inv[i] = ull(inv[i - 1]) * inv[i] % MOD;
    }
}
unsigned pow(unsigned a, unsigned b) {
    unsigned res = 1;
    while (b) {
        if (b & 1) {
            res = ull(res) * a % MOD;
        }
        a = ull(a) * a % MOD;
        b >>= 1;
    }
    return res;
}
main() {
    ios::sync_with_stdio(false);
    cin >> n >> k;
    init();
    for (unsigned i = 1; i <= k + 2; i++) {
        f[i] = (f[i - 1] + pow(i, k)) % MOD;
    }
    if (n <= k + 2) {
        cout << f[n];
        return 0;
    }
    unsigned t = 1;
    for (unsigned i = 1; i <= k + 2; i++) {
        t = ull(t) * (n - i) % MOD;
    }
    unsigned ans = 0;
    for (unsigned i = 1; i <= k + 2; i++) {
        ans = (ans + ((ll(ull(f[i]) * t % MOD * pow(n - i, MOD - 2) % MOD * inv[i - 1] % MOD * inv[k + 2 - i] % MOD) * (((k + 2 - i) & 1) ? -1 : 1)) + MOD) % MOD) % MOD;
    }
    cout << ans;
}