拉格朗日插值法

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

例题

The sum of the K-th Powers

给出两个非负整数 $n$n，$k$k，求
$\sum_{i=1}^ni^k$i=1∑n​ik
答案对 $10^9+7$109+7 取模。
$n\le10^9$n≤109，$k\le10^6$k≤106。

题解

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

代码

#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;
}