适用范围:给出n个点$(x_i,y_i)$,过这n个点能够确定一个最高n-1次的多项式$f(x)$,求$f(k)$

做法:如图所示,我们将每一个点$(x_i,y_i)$在x轴上的投影$(x_i,0)$记为$H_i$。对每一个i,我们选择一个点集$lbrace P_ibrace cup lbrace H_j vert 1 le ile n, j
eq ibrace$,求过这n个点的最高n-1次的多项式$g_i(x)$。这样我们就得到了n个$g_i(x)$,它们都在各自对应的$x_i$处的值为$y_i$,并且在其它$x_j(i
eq j)$处值为0

拉格朗日插值法学习笔记-编程知识网

很容易就能够构造出$g_i(x)$的表达式:

$$g_i(x)=y_i*prod_{j
eq i}frac{x-x_j}{x_i-x_j}$$

显然最后有:

$$f(x)=sum_{i=1}^{n}g_i(x)=sum _{i=1}^{n}y_i*prod_{j
eq i}frac{x-x_j}{x_i-x_j}$$

由于只用求$f(k)$的值,代入得$f(k)=sum _{i=1}^{n}y_i*prod_{j
eq i}frac{k-x_j}{x_i-x_j}$

例题:

luogu 4781 [模板]拉格朗日插值

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>

using namespace std;

typedef long long ll;

const ll mod = 998244353;
const int N = 2010;

int n;
ll k, x[N], y[N];

ll power(ll a, ll n)
{
    ll res = 1;
    while (n) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

int main()
{
    scanf("%d%lld", &n, &k);
    for (int i = 1; i <= n; i++) scanf("%lld%lld", &x[i], &y[i]);
    ll res = 0;
    for (int i = 1; i <= n; i++) {
        ll ta = y[i] % mod, tb = 1;
        for (int h = 1; h <= n; h++) {
            if (h == i) continue;
            ta = (ta * ((k - x[h]) % mod + mod) % mod) % mod;
            tb = (tb * ((x[i] - x[h]) % mod + mod) % mod) % mod;
        }
        res = (res + ta * power(tb, mod - 2) % mod) % mod;
    }
    printf("%lld
", res);
    return 0;
}

[模板]拉格朗日插值

Codeforces – 622F The Sum of the k-th Powers

题意:求$sum_{i=1}^{n}i^k,1leq nleq 10^9,0leq k leq 10^6$

思路:首先有一个结论,$sum_{i=1}^{n}i^k$为k+1阶多项式,所以我们只需要暴力算出$f(n)=sum_{i=1}^{n}i^k$的前k+2项,然后用拉格朗日插值法求第n项即可

下面简单证明一下$sum_{i=1}^{n}i^k$为k+1阶多项式

对于一个数列${a_n}$来说,把数列${a_n}$的元素两两做差得到数列${b_n}$,我们称数列${b_n}$为数列${a_n}$的一阶阶差数列,如果再将数列${b_n}$的元素两两做差得到数列${c_n}$,我们称数列${c_n}$为数列${a_n}$的一阶阶差数列,依次类推定义p阶阶差数列


定理:数列${a_n}$是一个p阶等差数列的充要条件是数列的通项$a_n$为n的一个p次多项式

证明:设$f(x)=sum_{i=0}^{n}u_i x^i$,令${b_n}$为${a_n}$的一阶阶差数列

那么$Delta f(x)=f(x+1)-f(x)=sum_{i=0}^{n}u_i (x+1)^i-sum_{i=0}^{n}u_i x^i$

我们只考虑$x^n$有$Delta f(x)=u_n(x+1)^n-u_n x^n$

将$u_n(x+1)^n$二项式展开,仅考虑$x^n$有$Delta f(x)=u_n x^n-u_n x^n=0$

所以每次差分后多项式的最高次数会减1,p次差分后,变为常数项,此时多项式的最高次数为0,定理成立

显然,如果一个多项式k次差分之后变为0,那么这个多项式的最高次数为k-1


我们将$sum_{i=1}^{n}i^k$差分后得$1^k,2^k,3^kdots n^k$

显然$1^k,2^k,3^kdots n^k$为k阶多项式,所以$sum_{i=1}^{n}i^k$为k+1阶多项式

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
 
using namespace std;
 
typedef long long ll;
 
const int N = 2000010;
const ll mod = 1000000007;
 
int n, k;
ll a[N], f[N], nf[N], now;
 
ll power(ll a, ll n)
{
    ll res = 1;
    while (n) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}
 
ll solve()
{
    now = f[0] = nf[0] = 1;
    for (int i = 1; i <= k + 2; i++) {
        now = now * (n - i) % mod;
        f[i] = f[i - 1] * i % mod;
        nf[i] = -nf[i - 1] * i % mod;
    }
    ll res = 0;
    for (int i = 1; i <= k + 2; i++) {
        ll t = power(f[i - 1] * nf[k + 2 - i] % mod, mod - 2);
        res = (res + a[i] * now % mod * power(n - i, mod - 2) % mod * t % mod) % mod;
        res = (res + mod) % mod;
    }
    return res;
}
 
int main()
{
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    scanf("%d%d", &n, &k);
    for (int i = 1; i <= min(k + 2, n); i++)
        a[i] = (a[i - 1] + power(i, k)) % mod;
    if (n <= k + 2) printf("%lld
", a[n]);
    else printf("%lld
", solve());
    return 0;
}

The Sum of the k-th Powers

luogu 4593 [TJOI2018]教科书般的亵渎

题意:求$sum_{i=0}^m f(n-a_i)-sum_{i=0}^{m-1} sum_{j=i+1}^{m}(a_j-a_i)^m+1$,其中$f(n)=sum_{i=1}^{n} i^{m+1}$

思路:后半部分$sum_{i=0}^{m-1} sum_{j=i+1}^{m}(a_j-a_i)^m+1$可以直接暴力求解

$f(n)=sum_{i=1}^{n} i^{m+1}$为m+2阶多项式,考虑到m不是很大,所以我们可以插入m+3个值后,暴力对每一个$a_i$求一次f(n-a_i),最后相加即可

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>

using namespace std;

typedef long long ll;

const int N = 60;
const ll mod = 1000000007;

int T;
ll n, m, now;
ll a[N], c[N], f[N], nf[N];

ll power(ll a, ll n)
{
    ll res = 1;
    while (n) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

ll solve(ll n)
{
    now = f[0] = nf[0] = 1;
    for (int i = 1; i <= m + 3; i++) {
        now = now * (n - i) % mod;
        f[i] = f[i - 1] * i % mod;
        nf[i] = -nf[i - 1] * i % mod;
    }
    ll res = 0;
    for (int i = 1; i <= m + 3; i++) {
        ll t = power(f[i - 1] * nf[m + 3 - i] % mod, mod - 2);
        res = (res + c[i] * now % mod * power(n - i, mod - 2) % mod * t % mod) % mod;
        res = (res + mod) % mod;
    }
    return res;
}

ll calc(ll x)
{
    if (x <= m + 3) return c[x];
    return solve(x);
}

int main()
{
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    scanf("%d", &T);
    while (T--) {
        scanf("%lld%lld", &n, &m);
        for (int i = 1; i <= m; i++) scanf("%lld", &a[i]);
        sort(a + 1, a + m + 1);
        for (int i = 1; i <= min(m + 3, n); i++)
            c[i] = (c[i - 1] + power(i, m + 1)) % mod;
        ll res = 0;
        for (int i = 0; i <= m; i++) res = (res + calc(n - a[i])) % mod;
        for (int i = 0; i <= m - 1; i++) {
            for (int h = i + 1; h <= m; h++) {
                ll t = power(a[h] - a[i], m + 1);
                res = ((res - t) % mod + mod) % mod;
            }
        }
        printf("%lld
", res);
    }
    return 0;
}

[TJOI2018]教科书般的亵渎

bzoj 3453 tyvj 1858 XLkxc

题意:求$sum_{i=0}^n sum_{j=1}^{a+i*d} sum_{x=1}^{j}x^k$

思路:$f(n)=sum_{x=1}^{n}x^k$,为k+1阶多项式

$g(n)=sum_{i=1}^{n} sum_{j=1}^{i} j^k$,g(n)差分后的结果为f(n),所以g(n)为k+2阶多项式

$res(n)=sum_{i=0}^{n}g(a+i*d)$,res(n)进行k+5次差分后为0,所以res(n)为k+4阶多项式

然后插值求解即可

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>

using namespace std;

typedef long long ll;

const int N = 200;
const ll mod = 1234567891;

int T, k;
ll a, n, d, now, fc[N], nfc[N];
ll f[N], g[N];

ll power(ll a, ll n)
{
    ll res = 1;
    while (n) {
        if (n & 1) res = res * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return res;
}

ll solve(ll n, int m, ll *c)
{
    now = fc[0] = nfc[0] = 1;
    for (int i = 1; i <= m; i++) {
        now = now * (n - i) % mod;
        fc[i] = fc[i - 1] * i % mod;
        nfc[i] = -nfc[i - 1] * i % mod;
    }
    ll res = 0;
    for (int i = 1; i <= m; i++) {
        ll t = power(fc[i - 1] * nfc[m - i] % mod, mod - 2);
        res = (res + c[i] * now % mod * power(n - i, mod - 2) % mod * t % mod) % mod;
        res = (res + mod) % mod;
    }
    return res;
}

ll calc(ll n, int m, ll *c)
{
    if (n <= m) return c[n];
    return solve(n, m, c);
}

int main()
{
    // freopen("in.txt", "r", stdin);
    // freopen("out.txt", "w", stdout);
    scanf("%d", &T);
    while (T--) {
        scanf("%d%lld%lld%lld", &k, &a, &n, &d);
        for (int i = 0; i <= k + 3; i++) f[i] = (f[i - 1] + power(i, k)) % mod;
        for (int i = 1; i <= k + 3; i++) f[i] = (f[i] + f[i - 1]) % mod;
        for (int i = 0; i <= k + 5; i++) g[i] = calc((a + i * d) % mod, k + 3, f);
        for (int i = 1; i <= k + 5; i++) g[i] = (g[i] + g[i - 1]) % mod;
        printf("%lld
", calc(n, k + 5, g));
    }
    return 0;
}

tyvj 1858 XLkxc

参考:https://oi-wiki.org/math/poly/lagrange/