拉格朗日插值

Lagrange插值

求拉格朗日插值需要注意的几个点，求f(x)之前先确定好需要几个点的坐标，求得所有的y[i]，这个只需要求一次，pre[]和suf[]在每求一次f(x)都需要计算一次，然后每一项y[i]需乘以pre[i-1]*suf[i+1]*inv(fac[i-1])*inv(fac[n+1-i])，每一项乘完之后要去mod，最后在乘以(-1)^{n+1-i}确定符号，最后加到答案里

inline ll f(ll x){
    ll ans = 0;
    pre[0] = suf[k+2+1] = 1;
    for (int i=1; i<=k+2; i++)  pre[i] = pre[i-1]*(x-i)%mod;
    for (int i=k+2; i>=1; i--)  suf[i] = suf[i+1]*(x-i)%mod;
    for (int i=1; i<=k+2; i++){
        int flag = 1;
        if ((k+2-i)&1)  flag = -1;
        ans = (ans+y[i]*pre[i-1]%mod*suf[i+1]%mod*inv[i-1]%mod*inv[k+2-i]%mod*flag)%mod;
    }
    return (ans+mod)%mod;
}

对于一个n次多项式，可以用系数表示法表示，这样需要n+1个系数；也可以用点值表示法表示，需要n+1个横坐标不同的点。也就是说，n+1个点可以唯一确定一个最高次不超过n次的多项式，可以用n+1个点来求得一个多项式。

给定n个点，确定一个多项式。
若要根据这n个点求得多项式的所有系数，则需要时间复杂度为O(n^3)的高斯消元法来解决，且存在精度问题。
利用拉格朗日插值法可以在O(n^2)的复杂度求得多项式在某个位置的值f(x).设已经给出了n+1个点(xi, yi).
公式：f(x) = \sum_{i=0}^n y[i]*(\prod_{j=0,j \neq i}^n\frac{x-x[j]}{x[i]-x[j]})

当这n+1个点的横坐标毫无规律时，只能O(n^2)求，对于每个y[i]，先把后面乘积符号的分子全都乘起来，分母也乘起来，然后分母去逆元即可。
当n+1个点的横坐标连续时，假设是[1, n+1]，则可以优化成O(n)来处理。
f(x) = \sum_{i=0}^ny[i]*(\prod_{j=0,j \neq i}^n\frac{x-x[j]}{x[i]-x[j]}) = \sum_{i=0}^ny[i]*(\prod_{j=0, j \neq i}^n\frac{x-j}{i-j}) =\sum_{i=0}^ny[i]*(\frac{pre[i-1]*suf[i+1]}{fac[i-1]*fac[n-i+1]*(-1)^{n-i+1}})

pre[i] = pre[i-1]*(x-i), pre[0] = 1
suf[i] = suf[i+1]*(x-i), suf[n+2] = 1

pre[]和suf[]都可以线性求得，阶乘fac[]提前维护即可。

求拉格朗日插值需要注意的几个点，求f(x)之前先确定好需要几个点的坐标，求得所有的y[i]，这个只需要求一次，pre[]和suf[]在每求一次f(x)都需要计算一次，然后每一项y[i]需乘以pre[i-1]*suf[i+1]*inv(fac[i-1])*inv(fac[n+1-i])，每一项乘完之后要去mod，最后在乘以(-1)^{n+1-i}确定符号，最后加到答案里

一个经典应用：求式子f(x) = \sum_{i=0}^x i^k.
这是一个k+1次多项式，证明略。所以求得前k+2个点后，可以O(k)求得任意一个f(x).当前面多一个求和符号时，最高次+1，多几个就+几。

P4781 拉格朗日插值

普通拉格朗日插值，O(n^2)求法。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int mod = 998244353, maxn = 2003;

ll n, k;
ll x[maxn], y[maxn];
inline void extgcd(ll a, ll b, ll &x, ll &y){
    if (b==0){
        x = 1, y = 0;
        return;
    }
    extgcd(b, a%b, y, x);
    y -= x*(a/b);
}
inline ll inv(ll a){
    ll x, y;
    a = (a%mod+mod)%mod;
    extgcd(a, mod, x, y);
    return (x%mod+mod)%mod;
}
int main(){
    ll ans = 0;
    cin >> n >> k;
    for (int i=1; i<=n; i++)    cin >> x[i] >> y[i];
    for (int i=1; i<=n; i++){
        for (int j=1; j<=n; j++)    if (j!=i)
            y[i] = y[i]*(k-x[j])%mod;
        ll tmp = 1;
        for (int j=1; j<=n; j++)    if (j!=i)
            tmp = tmp*(x[i]-x[j])%mod;
        ans += y[i]*inv(tmp)%mod;
    }
    ans = (ans%mod+mod)%mod;
    cout << ans;
    return 0;
}

CF622F The Sum of the k-th Powers

求式子: \sum_{i=0}^n i^k.

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const ll maxn = 1e6+5, mod = 1e9+7;
int n, k;
ll fac[maxn], pre[maxn], suf[maxn], y[maxn], inv[maxn];
ll pow_mod(ll a, ll b){
    a %= mod, b %= mod-1;
    if (!a) return 0;
    ll ans = 1;
    while (b){
        if (b&1)    ans = ans*a%mod;
        b >>= 1;
        a = a*a%mod;
    }
    return ans;
}
inline ll f(ll x){
    ll ans = 0;
    pre[0] = suf[k+2+1] = 1;
    for (int i=1; i<=k+2; i++)  pre[i] = pre[i-1]*(x-i)%mod;
    for (int i=k+2; i>=1; i--)  suf[i] = suf[i+1]*(x-i)%mod;
    for (int i=1; i<=k+2; i++){
        int flag = 1;
        if ((k+2-i)&1)  flag = -1;
        ans = (ans+y[i]*pre[i-1]%mod*suf[i+1]%mod*inv[i-1]%mod*inv[k+2-i]%mod*flag)%mod;
    }
    return (ans+mod)%mod;
}
int main(){
    cin >> n >> k;
    for (int i=1; i<=k+2; i++)  y[i] = (y[i-1]+pow_mod(i, k))%mod;
    fac[0] = 1, inv[0] = 1;
    for (int i=1; i<=k+2; i++)  fac[i] = fac[i-1]*i%mod, inv[i] = pow_mod(fac[i], mod-2);
    cout << f(n);
    return 0;
}

P4593 教科书般的亵渎

看题都看了半天。。懒得写了。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const ll mod = 1e9+7;
ll n, m, k;
ll a[55];
ll y[55], pre[55], suf[55], fac[55];
inline ll pow_mod(ll a, ll b){
    a %= mod;
    if (!a) return 0;
    b %= mod-1;
    ll ans = 1;
    while (b){
        if (b&1)    ans = ans*a%mod;
        b >>= 1;
        a = a*a%mod;
    }
    return ans;
}
inline ll inv(ll a){
    a = (a%mod+mod)%mod;
    return pow_mod(a, mod-2);
}
void init(){
    y[1] = 1;
    for (int i=2; i<=k+2; i++)  y[i] = (y[i-1]+pow_mod(i, k))%mod;
    return;
}
ll f(ll x){
    x %= mod;
    ll ans = 0;
    pre[0] = suf[k+2+1] = 1;
    for (int i=1; i<=k+2; i++)  pre[i] = pre[i-1]*(x-i)%mod;
    for (int i=k+2; i>=1; i--)  suf[i] = suf[i+1]*(x-i)%mod;
    for (int i=1; i<=k+2; i++){
        ans = (ans+y[i]*pre[i-1]%mod*suf[i+1]%mod*inv(fac[i-1])%mod*inv(fac[k+2-i])%mod*(((k+2-i)&1)?-1:1))%mod;
    }
    return (ans+mod)%mod;
}
int main(){
    int t;
    fac[0] = fac[1] = 1;
    for (int i=2; i<=53; i++)   fac[i] = fac[i-1]*i%mod;
    scanf ("%d", &t);
    while (t--){
        ll ans = 0;
        scanf ("%lld %lld", &n, &m);    k = m+1;
        for (int i=1; i<=m; i++)    scanf ("%lld", a+i);
        init();
        sort (a+1, a+m+1);
        for (int i=1; i<=m; i++){
            ans = (ans+f(n))%mod;
            for (int j=i; j<=m; j++)    ans = (ans-pow_mod(a[j], k))%mod;
            for (int j=i+1; j<=m; j++)    a[j] -= a[i];
            n -= a[i];
        }
        ans = (ans+f(n))%mod;
        printf ("%lld\n", (ans%mod+mod)%mod);
    }
    return 0;
}

LOJ6024 XLKxc

给定k, a, n, d, 求式子：\sum_{i=0}^n \sum_{j=1}^{a+id} \sum_{l=1}^{j} l^k (\%p)
对式子进行转化：\sum_{i=0}^n \sum_{j=1}^{a+id} \sum_{l=1}^{j} l^k (\%p) = \sum_{i=0}^n \sum_{j=1}^{a+id} f(j) (\%p) = \sum_{i=0}^{n}g(a+id) = h(n)
根据多项式的定义有：
f(x) = \sum_{i=1}^{x}i^k是k+1次多项式。
g(x) = \sum_{j=1}^{x}f(x)是k+2次多项式
h(x) = \sum_{i=0}^{x}g(x)是k+3次多项式
使用拉格朗日插值法可在O(k^2)复杂度内解决问题。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const ll mod = 1234567891;

ll k, a, n, d;
ll fac[150], suf[150], pre[150];
ll y[150], y2[150], y3[150];
inline ll pow_mod(ll a, ll b){
    a = (a%mod+mod)%mod;
    if (!a) return 0;
    b %= mod-1;
    ll ans = 1;
    while (b){
        if (b&1)    ans = ans*a%mod;
        b >>= 1;
        a = a*a%mod;
    }
    return ans;
}
inline ll inv(ll a){
    return pow_mod(a, mod-2);
}

inline ll f(ll x){
    ll ans = 0;
    x %= mod;
    pre[0] = suf[k+2+1] = 1;
    for (int i=1; i<=k+2; i++)  pre[i] = pre[i-1]*(x-i)%mod;
    for (int i=k+2; i>=1; i--)  suf[i] = suf[i+1]*(x-i)%mod;
    for (int i=1; i<=k+2; i++){
        int flag = 1;
        if ((k+2-i)&1)    flag = -1;
        ans = (ans+y[i]*pre[i-1]%mod*suf[i+1]%mod*inv(fac[i-1])%mod*inv(fac[k+2-i])%mod*flag)%mod;
    }
    return (ans%mod+mod)%mod;
}
inline ll g(ll x){
    ll ans = 0;
    x %= mod;
    pre[0] = suf[k+3+1] = 1;
    for (int i=1; i<=k+3; i++)  pre[i] = pre[i-1]*(x-i)%mod;
    for (int i=k+3; i>=1; i--)  suf[i] = suf[i+1]*(x-i)%mod;
    for (int i=1; i<=k+3; i++){
        int flag = 1;
        if ((k+3-i)&1)  flag = -1;
        ans = (ans+y2[i]*pre[i-1]%mod*suf[i+1]%mod*inv(fac[i-1])%mod*inv(fac[k+3-i])%mod*flag)%mod;
    }
    return (ans%mod+mod)%mod;
}
inline ll h(ll x){
    ll ans = 0;
    x %= mod;
    pre[0] = suf[k+4+1] = 1;
    for (int i=1; i<=k+4; i++)  pre[i] = pre[i-1]*(x-i)%mod;
    for (int i=k+4; i>=1; i--)  suf[i] = suf[i+1]*(x-i)%mod;
    for (int i=1; i<=k+4; i++){
        int flag = 1;
        if ((k+4-i)&1)  flag = -1;
        ans = (ans+y3[i]*pre[i-1]%mod*suf[i+1]%mod*inv(fac[i-1])%mod*inv(fac[k+4-i])%mod*flag)%mod;
    }
    return (ans%mod+mod)%mod;
}
int main(){
    int t;
    fac[0] = 1;
    for (int i=1; i<=130; i++)  fac[i] = fac[i-1]*i%mod;
    scanf ("%d", &t);
    while (t--){
        scanf ("%lld %lld %lld %lld", &k, &a, &n, &d);
        y[0] = 0, y2[0] = 0, y3[0] = 0;
        for (int i=1; i<=k+2; i++)  y[i] = (y[i-1]+pow_mod(i, k))%mod;           //y只要求1次即可，之后会重复用到
        for (int i=1; i<=k+3; i++)  y2[i] = (y2[i-1]+f(i))%mod;
        y3[0] = g(a);
        for (int i=1; i<=k+4; i++)  y3[i] = (y3[i-1]+g(a+i*d))%mod;
        printf ("%lld\n", h(n));
    }
    return 0;
}

ps: 对于一个多项式f(x)，前面加一个求和符号后，还是一个多项式，最高次+1.