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

【算法讲14：拉格朗日插值】

• 前言
• 引入
• ⌈ \lceil ⌈拉格朗日插值 ⌋ \rfloor ⌋
• 代码
• ⌈ \lceil ⌈拉格朗日插值 ⌋ \rfloor ⌋差分法运用
• 代码

前言

• 拉格朗日插值可以出的很难，对于一个 n n n 阶多项式，可以优化成 O ( n log ⁡ n ) O(n\log n) O(nlogn)

  但是需要 F F T FFT FFT 等知识点，博主比较菜就只能搞搞入门的了 (哭)

• 基础的高斯消元法可以 O ( n 3 ) O(n^3) O(n3) 求解 n n n 阶多项式，但是貌似不够用呀

  这时就出现了 O ( n 2 ) O(n^2) O(n2) 的拉格朗日插值！

  对于能求出 f ( i ) , ∀ i ∈ [ 1 , n ] f(i),\forall i\in[1,n] f(i),∀i∈[1,n]，我们用 差分法 还可以化简到 O ( n ) O(n) O(n)

引入

• 【题目描述】【模板】拉格朗日插值 | 洛谷 P4781

  有 n n n 个不同点，给定他们的坐标 ( x i , y i ) (x_i,y_i) (xi​,yi​)

  由于他们的横坐标都不相同，因此他们可以唯一确定一个 n − 1 n-1 n−1 阶的多项式 P ( x ) P(x) P(x)

  给定 k k k ，你的任务是求出 P ( k ) P(k) P(k) 取模 998244353 998244353 998244353

• 【数据范围】

  1 ≤ n ≤ 2 × 1 0 3 1\le n\le 2\times10^3 1≤n≤2×103

  1 ≤ x i , y i , k < 998244353 1\le x_i,y_i,k<998244353 1≤xi​,yi​,k<998244353

⌈ \lceil ⌈拉格朗日插值 ⌋ \rfloor ⌋

• 我们知道多项式可以表示为 P ( x ) = ∑ a i x i P(x)=\sum a_ix^i P(x)=∑ai​xi。

  很神奇的，我们想构造出一个函数来表示这个 P ( x ) P(x) P(x)

• 我们设一个函数：

  l i ( x ) = ∏ j = 0 , j ≠ i n x − x j x i − x j l_i(x)=\prod_{j=0,j\ne i}^n\cfrac{x-x_j}{x_i-x_j} li​(x)=j=0,j​=i∏n​xi​−xj​x−xj​​

  这个函数有什么性质呢？有一下两个性质：

  （1） ∀ i ∈ [ 0 , n ] ， l i ( x i ) = 1 \forall i\in[0,n]，l_i(x_i)=1 ∀i∈[0,n]，li​(xi​)=1

  （2） ∀ i , j ∈ [ 0 , n ] , i ≠ j ， l i ( x j ) = 0 \forall i,j\in[0,n],i\ne j，l_i(x_j)=0 ∀i,j∈[0,n],i​=j，li​(xj​)=0

  接下来，我们就能直接构造出 P ( x ) P(x) P(x)了：

  P ( x ) = ∑ i = 0 n y i l i ( x ) P(x)=\sum_{i=0}^n y_il_i(x) P(x)=i=0∑n​yi​li​(x)

  为啥这个多项式保证是正确的呢？因为我们有：

  ∀ i ∈ [ 0 , n ] ， P ( x i ) = y i \forall i\in[0,n]，P(x_i)=y_i ∀i∈[0,n]，P(xi​)=yi​，保证成立。

代码

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

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 n、k n、k，求：

  ∑ i = 1 n i k \sum_{i=1}^ni^k i=1∑n​ik

• 【数据范围】

  1 ≤ n ≤ 1 0 9 1\le n\le 10^9 1≤n≤109

  1 ≤ k ≤ 1 0 6 1\le k\le 10^6 1≤k≤106

• 【思路】

  这下阶为 1 e 6 1e6 1e6了， O ( k 2 ) O(k^2) O(k2) 已经不适用了。

  我们如果算出前面一些 x i x_i xi​ 和 y i y_i yi​，想快速算 l i ( n ) l_i(n) li​(n)，怎么办呢？

  如果我们令 x i = i x_i=i xi​=i，那么式子会变成下面的样子：

  l i ( n ) = ∏ j = 1 , j ≠ i k n − j i − j l_i(n)=\prod_{j=1,j\ne i}^k\cfrac{n-j}{i-j} li​(n)=j=1,j​=i∏k​i−jn−j​

  （1）对于分子，我们只要预处理出一个前缀积和后缀积，即可快速算。

  令 p r e ( i ) = ∏ j = 1 i n − j pre(i)=\prod_{j=1}^i n-j pre(i)=∏j=1i​n−j， s u f ( i ) = ∏ j = i k n − j suf(i)=\prod_{j=i}^k n-j suf(i)=∏j=ik​n−j

  那么分子就是 p r e ( i − 1 ) × s u f ( i + 1 ) pre(i-1)\times suf(i+1) pre(i−1)×suf(i+1)

  （2）对于分母，容易想到就是一个 断开的阶乘，就是 f a c [ i − 1 ] × f a c [ k − i ] fac[i-1]\times fac[k-i] fac[i−1]×fac[k−i]

  不过要注意一下，因为当 j > i j>i j>i 时，后面每一项都为负数。若 k − i k-i k−i是一个奇数，那么后面会多一个负号。

  P ( n ) = ∑ i = 1 k y i × p r e ( i − 1 ) × s u f ( i + 1 ) f a c [ i ] × f a c [ k − i ] × ( − 1 ) k − i P(n)=\sum_{i=1}^ky_i\times\cfrac{pre(i-1)\times suf(i+1)}{fac[i]\times fac[k-i]}\times(-1)^{k-i} P(n)=i=1∑k​yi​×fac[i]×fac[k−i]pre(i−1)×suf(i+1)​×(−1)k−i

  这下，时间复杂度变为线性的了！

代码

• 342 M s / 2000 M s 342Ms/2000Ms 342Ms/2000Ms

  时间复杂度： O ( k ) O(k) 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;
}