拉格朗日插值合辑

什么是插值？

在最开始，我们知道，两个点能确定一个一元一次函数，三个点能确定一个一元二次函数…由此推广，n+1个点能确定一个一元n次函数。插值法就是利用f(x)在某区间插入若干点的函数值，用这些点的值描绘出一个特定的近似函数曲线，用这条函数曲线能表示出其他值。如果这个函数是一个多项式，那么就称它为插值多项式。

拉格朗日插值有啥用？

拉格朗日插值就是一个n次多项式的插值。
n+1个点能确定一个最高次为n次的多项式，题目中通常的将此当做一个辅助工具，通常是在处理出一个函数以后，一般这个函数能够确定最高次数，题目直接或者间接的给出了n+1个点，此时需要知道一个自变量很大的函数值，拉格朗日插值就派上用场了。

拉格朗日插值

关键是围绕这个公式展开的：$$f(x_i)=\sum _{i=0}^n{y}_i\prod _{j=0,i\ne j}^n\frac{x_i-x[j]}{x[i]-x[j]}$$
从这个公式中也能看出，n+1个点能近似表示一个n次函数，并且对于任意一个点$x_i$，能通过这n+1个点得到$f(x_i)$，这就是它的灵魂所在，这个算法能看出来复杂度是$O(n^2)$的

在给出的n+1个点x值是连续的情况

在这个情况下，我们能对公式化简，使其到达$O(n)$的复杂度。$$f(x_i)=\sum _{i=0}^n{y}_i\prod _{j=0,i\ne j}^n\frac{x_i-j}{i-j}$$
对于分子而言处理处它的前缀积和后缀积。
前缀积：$pr{e}_i={\prod }_{j=0}^i{x}_i-j$，后缀积：$su{f}_i={\prod }_{j=i}^n{x}_i-j$
分母的话能化成两个阶乘相乘的形式，处理出阶乘逆元即可，于是就有：$$f(x_i)=\sum _{i=0}^n{y}_i\prod _{j=0,i\ne j}^n\frac{pr{e}_i\ast su{f}_i}{frac[i]\ast frac[n-i]}$$
分母中要注意：如果i是从0~n，那么分母是$frac[i]\ast frac[n-i]$，如果i是从1~n+1的话，这时分母记得要换成$frac[i-1]\ast frac[n-i]$。

一些例题

CF622F The Sum of the k-th Powers

题意：给出n,k，求${\sum }_{i=1}^n{i}^k。(1\le n\le 10^9,0\le k\le 10^6).$
做法：n为10^9，不能直接循环得出，若干k次项的前n项和为k+1次项式，其函数构造的点可以用$(1,1^k),(2,1^k+2^k),(3,1^k+2^k+3^k),...,(k+2,1^k+2^k+...+(k+2{)}^k)$这k+2个点，然后算f(n)就是我们要的答案。
代码：

#include  <algorithm>
#include  <iostream>
#include  <cstring>
#include  <vector>
#include  <cstdio>
#include  <queue>
#include  <cmath>
#include  <set>
#include  <map>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define _for(n,m,i) for (register int i = (n); i < (m); ++i)
#define _rep(n,m,i) for (register int i = (n); i <= (m); ++i)
#define lson rt<<1, l, mid
#define rson rt<<1|1, mid+1, r
#define PI acos(-1)
#define eps 1e-8
#define rint register int
#define F(x) ((x)/3+((x)%3==1?0:tb))
#define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)
using namespace std;
typedef long long LL;
typedef pair<LL, int> pli;
typedef pair<int, int> pii;
typedef pair<double, int> pdi;
typedef pair<LL, LL> pll;
typedef pair<double, double> pdd;
typedef map<int, int> mii;
typedef map<char, int> mci;
typedef map<string, int> msi;
template<class T>
void read(T &res) {
  int f = 1; res = 0;
  char c = getchar();
  while(c < '0' || c > '9') { if(c == '-') f = -1; c = getchar(); }
  while(c >= '0' && c <= '9') { res = res * 10 + c - '0'; c = getchar(); }
  res *= f;
}
const int ne[8][2] = {1, 0, -1, 0, 0, 1, 0, -1, -1, -1, -1, 1, 1, -1, 1, 1};
const int INF = 0x3f3f3f3f;
const int N = 1e6+10;
const LL Mod = 1e9+7;
const int M = 1e6+10;
LL ksm(LL x, LL n) {
  LL res = 1;
  while(n) {
    if(n & 1) res = res * x % Mod;
    x = x * x % Mod;
    n >>= 1;
  }
  return res;
}
LL fac[N+1], inv[N+1];
void init() {
  fac[0] = 1;
  for(int i = 1; i <= N; i++) {
    fac[i] = fac[i-1] * i % Mod;
  }
  inv[N] = ksm(fac[N], Mod-2);
  for(int i = N-1; i >= 0; i--) {
    inv[i] = (inv[i+1] * (i + 1)) % Mod;
  }
}
LL sumy[N], pre[N], suf[N];
int main() {
  LL n, k; init();
  scanf("%lld %lld", &n, &k);

  pre[0] = suf[k+3] = 1;
  _rep(1, k+2, i) pre[i] = pre[i-1] * (n-i) % Mod;
  for(int i = k+2; i >= 1; --i) suf[i] = suf[i+1] * (n-i) % Mod;

  LL ans = 0, y = 0;
  _rep(1, k+2, i) {
    y = (y + ksm(i, k)) % Mod;
    ans = (ans + y*pre[i-1]%Mod*suf[i+1]%Mod*inv[i-1]%Mod*inv[k+2-i]%Mod*(((k+2-i)&1)?-1:1)) % Mod;
  }
  printf("%lld\n", (ans+Mod)%Mod);
  return 0;
}

bzoj3453 tyvj 1858 XLkxc

题意：
给定 $k,a,n,d,p,$

$f(i)=1^k+2^k+3^k+...+i^k，$

$g(x)=f(1)+f(2)+f(3)+...+f(x)，$

求$F(x)=(g(a)+g(a+d)+g(a+2d)+...+g(a+nd))modp$

做法：首先$$g(x)=\sum _{j=1}^x\sum _{l=1}^j{l}^k$$

$$F(x)=\sum _{i=0}^n\sum _{j=1}^{a+i\ast d}\sum _{l=1}^j{l}^k$$
f函数是一个k+1项式函数，那么作为f函数的前缀和的g函数则是一个k+2项式函数，先用前缀和处理出g数组，得到k+3个点，再分别用a,a+d,…,a+n*d作为$x_i$代入，用插值求出$g(x_i)$，用插值出来的函数做一次前缀和再来用n插值，就得到了结果，就是插值套插值。

代码：

#include<bits/stdc++.h>
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
#define _for(n,m,i) for (register int i = (n); i < (m); ++i)
#define _rep(n,m,i) for (register int i = (n); i <= (m); ++i)
#define lson rt<<1, l, mid
#define rson rt<<1|1, mid+1, r
#define PI acos(-1)
#define eps 1e-8
#define rint register int
using namespace std;
typedef long long LL;
typedef pair<LL, int> pli;
typedef pair<int, int> pii;
typedef pair<double, int> pdi;
typedef pair<LL, LL> pll;
typedef pair<double, double> pdd;
typedef map<int, int> mii;
typedef map<char, int> mci;
typedef map<string, int> msi;
template<class T>
void read(T &res) {
  int f = 1; res = 0;
  char c = getchar();
  while(c < '0' || c > '9') { if(c == '-') f = -1; c = getchar(); }
  while(c >= '0' && c <= '9') { res = res * 10 + c - '0'; c = getchar(); }
  res *= f;
}
const int ne[8][2] = {1, 0, -1, 0, 0, 1, 0, -1, -1, -1, -1, 1, 1, -1, 1, 1};
const int INF = 0x3f3f3f3f;
const LL N = 210;
const LL Mod = 1234567891;
const int M = 1e6+10;
LL ksm(LL x, LL n) {
  LL res = 1;
  while(n) {
    if(n & 1) res = res * x % Mod;
    x = x * x % Mod;
    n >>= 1;
  }
  return res;
}
LL fac[N+5], inv[N+5];
void init() {
  fac[0] = 1;
  for(int i = 1; i <= N; i++) {
    fac[i] = fac[i-1] * i % Mod;
  }
  inv[N] = ksm(fac[N], Mod-2);
  for(int i = N-1; i >= 0; i--) {
    inv[i] = (inv[i+1] * (i + 1)) % Mod;
  }
}
/*x取值连续，O(n)求法。
有n项，输入y数组，计算f[xi]是多少。*/
LL pre[N], suf[N];
LL Lagrange(LL n, LL *y, LL xi) {
  pre[0] = xi; suf[n+1] = 1;
  for(int i = 1; i <= n; ++i) pre[i] = pre[i-1]*(xi-i)%Mod;
  for(int i = n; i >= 0; --i) suf[i] = suf[i+1]*(xi-i)%Mod;
  LL ans = 0;
  _rep(0, n, i) {
    ans = (ans + y[i] * (i?pre[i-1]:1) % Mod * suf[i+1] % Mod * inv[i] % Mod * inv[n-i] % Mod * (((n-i)&1) ? -1: 1)) % Mod;
  }
  return (ans + Mod) % Mod;
}

LL f[N], g[N];
int main() {
  int T; scanf("%d", &T);
  LL k, a, n, d;
  init();
  while(T--) {
    scanf("%lld %lld %lld %lld", &k, &a, &n, &d);
    _rep(1, k+3, i) f[i] = (f[i-1] + ksm(i, k)) % Mod;
    _rep(1, k+3, i) f[i] = (f[i-1] + f[i]) % Mod;
    g[0] = Lagrange(k+2, f, a);
    _rep(1, k+3, i) {
      g[i] = (g[i-1] + Lagrange(k+2, f, (a+i*d) % Mod)) % Mod;
    }
    printf("%lld\n", Lagrange(k+3, g, n));
  }
  return 0;
}

最后推荐一篇优质的拉格朗日插值博客