【算法讲14：拉格朗日插值】拉格朗日插值入门 与 拉格朗日插值差分法

前言

• 拉格朗日插值可以出的很难，对于一个 $n$ 阶多项式，可以优化成 $O(n\log n)$
  但是需要 $FFT$ 等知识点，博主比较菜就只能搞搞入门的了 (哭)
• 基础的高斯消元法可以 $O(n^3)$ 求解 $n$ 阶多项式，但是貌似不够用呀
  这时就出现了 $O(n^2)$ 的拉格朗日插值！
  对于能求出 $f(i),\forall i\in [1,n]$，我们用 差分法 还可以化简到 $O(n)$

引入

• 【题目描述】【模板】拉格朗日插值 | 洛谷 P4781
  有 $n$ 个不同点，给定他们的坐标 $(x_i,y_i)$
  由于他们的横坐标都不相同，因此他们可以唯一确定一个 $n-1$ 阶的多项式 $P(x)$
  给定 $k$ ，你的任务是求出 $P(k)$ 取模 $998244353$
• 【数据范围】
  $1\le n\le 2\times 10^3$
  $1\le x_i,y_i,k<998244353$

$\lceil$拉格朗日插值$\rfloor$

• 我们知道多项式可以表示为 $P(x)=\sum a_i{x}^i$。
  很神奇的，我们想构造出一个函数来表示这个 $P(x)$
• 我们设一个函数：
  $$l_i(x)=\prod _{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}$$
  这个函数有什么性质呢？有一下两个性质：
  （1）$\forall i\in [0,n]，l_i(x_i)=1$
  （2）$\forall i,j\in [0,n],i\ne j，l_i(x_j)=0$
  接下来，我们就能直接构造出 $P(x)$了：
  $$P(x)=\sum _{i=0}^n{y}_i{l}_i(x)$$
  为啥这个多项式保证是正确的呢？因为我们有：
  $\forall i\in [0,n]，P(x_i)=y_i$，保证成立。

代码

• 时间复杂度：$O(n^2)$

const int MAX = 2e3+50;
const ll  MOD = 998244353;

ll qpow(ll a,ll n){/* */ll res = 1LL;while(n){if(n&1)res=res*a%MOD;a=a*a%MOD;n>>=1;}return res;}
ll inv(ll a){/* */return qpow(a,MOD-2);}

ll X[MAX],Y[MAX];

int main()
{
    int n;ll k;
    cin >> n >> k;
    for(int i = 0;i < n;++i){
        cin >> X[i] >> Y[i];
    }
    ll res = 0;
    for(int i = 0;i < n;++i){
        ll num = Y[i];
        ll den = 1;
        for(int j = 0;j < n;++j){
            if(i != j){
                num = num * (k - X[j]) % MOD;
                den = den * (X[i] - X[j]) % MOD;
            }
        }
        num = (num + MOD) % MOD;
        den = (den + MOD) % MOD;

        res = (res + num % MOD * inv(den) % MOD) % MOD;
    }
    cout << res;
    return 0;
}

$\lceil$拉格朗日插值$\rfloor$差分法运用

• 【题目描述】The Sum of the k-th Powers | CF 622F
  给定 $n、k$，求：
  $$\sum _{i=1}^n{i}^k$$
• 【数据范围】
  $1\le n\le 10^9$
  $1\le k\le 10^6$
• 【思路】
  这下阶为 $1e6$了，$O(k^2)$ 已经不适用了。
  我们如果算出前面一些 $x_i$ 和 $y_i$，想快速算 $l_i(n)$，怎么办呢？
  如果我们令 $x_i=i$，那么式子会变成下面的样子：
  $$l_i(n)=\prod _{j=1,j\ne i}^k\frac{n-j}{i-j}$$
  （1）对于分子，我们只要预处理出一个前缀积和后缀积，即可快速算。
  令 $pre(i)={\prod }_{j=1}^{i}n-j$，$suf(i)={\prod }_{j=i}^{k}n-j$
  那么分子就是 $pre(i-1)\times suf(i+1)$
  （2）对于分母，容易想到就是一个 断开的阶乘，就是 $fac[i-1]\times fac[k-i]$
  不过要注意一下，因为当 $j>i$ 时，后面每一项都为负数。若 $k-i$是一个奇数，那么后面会多一个负号。
  $$P(n)=\sum _{i=1}^k{y}_i\times \frac{pre(i-1)\times suf(i+1)}{fac[i]\times fac[k-i]}\times (-1{)}^{k-i}$$
  这下，时间复杂度变为线性的了！

代码

• $342Ms/2000Ms$
  时间复杂度：$O(k)$

/*
 _            __   __          _          _
| |           \ \ / /         | |        (_)
| |__  _   _   \ V /__ _ _ __ | |     ___ _
| '_ \| | | |   \ // _` | '_ \| |    / _ \ |
| |_) | |_| |   | | (_| | | | | |___|  __/ |
|_.__/ \__, |   \_/\__,_|_| |_\_____/\___|_|
        __/ |
       |___/
*/
const int MAX = 1e6+50;
const ll  MOD = 1e9+7;

ll qpow(ll a,ll n){/* */ll res = 1LL;while(n){if(n&1)res=res*a%MOD;a=a*a%MOD;n>>=1;}return res;}
ll inv(ll a){/* */return qpow(a,MOD-2);}


ll Y[MAX];
ll fac[MAX];
ll ivfac[MAX];
ll pre[MAX],suf[MAX];
void init(ll n,ll k){
    fac[0] = 1;
    pre[0] = 1;
    for(int i = 1;i <= k;++i){
        fac[i] = (fac[i-1] * i) % MOD;
        pre[i] = pre[i - 1] * (n - i) % MOD;
    }
    ivfac[k] = inv(fac[k]);
    suf[k+1] = 1;			/// 注意一下，不存在的话设为 1
    suf[k] = n - k;
    for(int i = k - 1;i >= 0;--i){
        ivfac[i] = ivfac[i + 1] * (i + 1) % MOD;
        suf[i] = suf[i + 1] * (n - i) % MOD;
    }
}
int main()
{
    ll n,k;
    cin >> n >> k;
    k += 2;			/// 因为是 k+2 阶的多项式
    init(n,k);
    ll sum = 0;
    for(int i = 1;i <= k;++i){
        sum  = (sum + qpow(i,k - 2)) % MOD;
        Y[i] = sum;
        if(n == i){
            cout << sum;
            return 0;
        }
    }
    ll res = 0;
    for(int i = 1;i <= k;++i){
        ll num = Y[i] * pre[i - 1] % MOD * suf[i + 1] % MOD;
        ll den = ivfac[i - 1] * ivfac[k - i] % MOD;
        if((k - i) % 2 == 1)den = -den;
        res = (res + num * den) % MOD;
    }
    res = (res + MOD) % MOD;
    cout << res;
    return 0;
}