【学习笔记】拉格朗日插值

整理的算法模板合集： ACM模板

P4781 【模板】拉格朗日插值

我们按照上面的思路，带入n个点求一下f(k)f(k)f(k)即可。时间复杂度O(n2)O(n^2)O(n2)。

注意本题还要求逆元，为了防止求逆元的时间复杂度影响整体的时间复杂度，所以我们分别计算出分子和分母，再将分子乘进分母的逆元，累加进最后的答案，时间复杂度的瓶颈就不会在求逆元上，总体的时间复杂度为O(n2)O(n^2)O(n2)。

代码实现：

本题需要累乘，会爆int，记得开long long

Code

// Problem: P4781 【模板】拉格朗日插值
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P4781
// Memory Limit: 125 MB
// Time Limit: 1000 ms
//
// Powered by CP Editor (https://cpeditor.org)#include <cstdio>
#include <iostream>
#include <cstring>
#include <algorithm>using namespace std;const int N = 500007;
typedef long long ll;
const int mod = 998244353;
ll n, m, k;struct Point
{ll x, y;
}A[N];ll qpow(ll a, ll b, ll c)
{ll res = 1;while(b) {if(b & 1) res = res * a % c;a = a * a % c;b >>= 1;}return res;
}ll inv(ll x) {return qpow(x, mod - 2, mod);}int main()
{scanf("%lld%lld", &n, &k);for(int i = 1; i <= n; ++ i) {scanf("%lld%lld", &A[i].x, &A[i].y);}ll ans = 0;for(int i = 1; i <= n; ++ i) {ll s1 = A[i].y % mod;ll s2 = 1ll;for(int j = 1; j <= n; ++ j) {if(i != j) {s1 = s1 * (k - A[j].x) % mod;s2 = s2 * (A[i].x - A[j].x) % mod;}}ans += s1 * inv(s2) % mod;}printf("%lld\n", (ans % mod + mod) % mod);return 0;
}

我们使用拉格朗日插值公式，对于一个数 kkk ，我们很容易在 O(n2)O(n^2)O(n2) 时间求得 F(k)F(k)F(k) 的数值，如果 xix_ixi 是连续的，我们甚至可以利用预处理在 O(n)O(n)O(n) 时间内得到 F(k)F(k)F(k) 的数值。
但是如果 xix_ixi 不连续，又有多组查询，就需要得到这个多项式的系数以保证求一个函数值的时间为 O(n)O(n)O(n) 。

重心拉格朗日插值法

我们对拉格朗日插值公式进行优化：

f(x)=∑i=1nyi∏j≠ix−xjxi−xjf(x)=\sum_{i=1}^n y_i\prod_{j\neq i}\dfrac {x-x_j} {x_i-x_j}f(x)=i=1∑nyij=i∏xi−xjx−xj

设 h=∏i=1nx−xih=\prod_{i=1}^n x-x_ih=i=1∏nx−xi
带入得：
=h∑i=1n∏j≠iyi(xi−xj)(x−xi)=h\sum_{i=1}^n \prod_{j\neq i}\dfrac {y_i} {(x_i-x_j)(x-x_i)}=hi=1∑nj=i∏(xi−xj)(x−xi)yi

设
ti=∏j≠iyixi−xjt_i=\prod_{j\neq i} \dfrac {y_i} {x_i-x_j}ti=j=i∏xi−xjyi

带入得：

=h∑i=1ntix−xi=h\sum_{i=1}^n \dfrac {t_i} {x-x_i}=hi=1∑nx−xiti

每次加入一个新点时，计算出它的 tit_iti ，并且更新别的点的 tit_iti ，时间复杂度 O(n)O(n)O(n)，加入 nnn 个点就是 4n^2$ 。

求解时， O(n)O(n)O(n) 求 hhh ，O(n)O(n)O(n) 求出那个 ∑\sum∑，时间复杂度 O(n)O(n)O(n)。

拉格朗日插值法求系数

首先观察式子，发现对于每个 iii , 上面部分是 (x−x1)×(x−x2)…×(x−xn)(x−xi)\cfrac{(x-x_1) \times (x-x_2)… \times (x-x_n)}{(x-x_i)}(x−xi)(x−x1)×(x−x2)…×(x−xn), 下面那部分和 yiy_iyi 都是常数。

所以可以 O(n2)O(n^2)O(n2) 处理出(x−x1)×(x−x2)…×(x−xn)(x-x_1) \times (x-x_2)… \times (x-x_n)(x−x1)×(x−x2)…×(x−xn) 这个 nnn 次多项式，然后通过模拟长除法，O(n)O(n)O(n)时间内可以得到(x−x1)×(x−x2)…×(x−xn)(x−xi)\cfrac{(x-x_1) \times (x-x_2)… \times (x-x_n)}{(x-x_i)}(x−xi)(x−x1)×(x−x2)…×(x−xn)，然后常系数直接 O(n)O(n)O(n) 暴力算出来，就得到了 xix_ixi 对应的多项式，最后把所有多项式加起来就得到了最终的系数。

// Problem: P4781 【模板】拉格朗日插值
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P4781
// Memory Limit: 125 MB
// Time Limit: 1000 ms
//
// Powered by CP Editor (https://cpeditor.org)#include <bits/stdc++.h>
using namespace std;
#define int long long
#define ll long long
const int N = 5007, mod = 998244353;int n, m, k;
int a[N], b[N], c[N], f[N], tmp[N], x[N], y[N];int qpow(int a, int b)
{int res = 1;while(b) {if(b & 1) res = 1ll * res * a % mod;a = 1ll * a * a % mod;b >>= 1;}return res % mod;
}int inv(int x)
{return qpow(x, mod - 2);
}void mul(ll *f, int len, ll t, int n) { //len为多项式的次数+1，f -> f*(x+t)for (int i = len; i > 0; -- i)tmp[i] = f[i], f[i] = f[i - 1];tmp[0] = f[0], f[0] = 0;for (int i = 0; i <= len; ++ i)f[i] = (f[i] + t * tmp[i]) % mod;
}
void dev(ll *f, ll *r, ll t, int n) { //f是被除多项式的系数，r = f/(x+t)for (int i = 0; i <= n; ++ i)tmp[i] = f[i];for (int i = n; i > 0; -- i) {r[i - 1] = tmp[i];tmp[i - 1] = (tmp[i - 1] - t * tmp[i] + mod) % mod;}return;
}//计算n次多项式 f_x = ∑{0~n-1} (a_i * x^i) 的 a_0 ~ a_{n - 1} (n次多项式系数)
void lagrange(int n) {  memset(a, 0, sizeof a);b[1] = 1, b[0] = -x[1];for (int i = 2; i <= n; ++i) {mul(b, i, -x[i], n);}//预处理(x-x1)*(x-x2)...*(x-xn)for (int i = 1; i <= n; ++i) {ll fz = 1;for (int j = 1; j <= n; ++j) {if (j == i)continue;fz = (fz * (x[i] - x[j]) + mod) % mod;}fz = qpow(fz, mod - 2);fz = fz * y[i] % mod; //得到多项式系数dev(b, c, -x[i], n);//得到多项式，保存在b数组for (int j = 0; j < n; ++j)a[j] = (a[j] + fz * c[j] % mod + mod) % mod;}
} signed main()
{scanf("%lld%lld", &n, &k);for(int i = 1; i <= n; ++ i) {scanf("%lld%lld", &x[i], &y[i]);}lagrange(n);int ans = 0;int kpow = 1;for(int i = 0; i < n; ++ i) {ans = (ans + kpow * a[i] % mod) % mod;kpow = kpow * k % mod;} printf("%lld\n", ans);return 0;
}

自然数k次幂的和

求 ∑i=1nik\sum_{i=1}^n i^k∑i=1nik ，n≤109,k≤106n \leq 10^9,k \leq 10^6n≤109,k≤106

Solution

n=k+2n=k+2n=k+2

// Problem: B. Morning Jogging
// Contest: Codeforces - Contest 2050 and Codeforces Round #718 (Div. 1 + Div. 2)
// URL: https://codeforces.com/contest/1517/problem/B
// Memory Limit: 256 MB
// Time Limit: 1000 ms
//
// Powered by CP Editor (https://cpeditor.org)#include <bits/stdc++.h>
using namespace std;const int N = 1100007, mod = 1e9 + 7;int read(){int s = 0, ne = 1; char c = getchar();while(c < '0' || c > '9') {if(c == '-') ne = -1; c = getchar();}while(c >= '0' && c <= '9') s = (s << 1) + (s << 3) + c - '0', c = getchar();return s * ne;
}int qpow(int a, int b)
{int res = 1;while(b) {if(b & 1) res = 1ll * res * a % mod;a = 1ll * a * a % mod; b >>= 1;} return res;
}int n, k;
int s[N], pre[N], suf[N], infact[N], fact[N], ans;void init()
{infact[0] = fact[0] = 1;for(int i = 1; i <= k + 10; ++ i)fact[i] = 1ll * fact[i - 1] * i % mod;infact[k + 10] = qpow(fact[k + 10], mod - 2);for(int i = k + 9; i >= 0; -- i)  infact[i] = 1ll * infact[i + 1] * (i + 1) % mod;infact[0] = 1;for(int i = 1; i <= k + 2; ++ i)s[i] = (s[i - 1] + qpow(i, k)) % mod;
}int main()
{scanf("%d%d", &n, &k);init();if(n <= k + 2) {printf("%d", s[n]);return 0;}pre[0] = 1;for(int i = 1; i <= k + 2; ++ i)pre[i] = 1ll * pre[i - 1] * (n - i) % mod;suf[k + 3] = 1;for(int i = k + 2; i; -- i)suf[i] = 1ll * suf[i + 1] * (n - i) % mod;for(int i = 1; i <= k + 2; ++ i) {s[i] = 1ll * s[i] * pre[i - 1] % mod * suf[i + 1] % mod * infact[i - 1] % mod * infact[k + 2 - i] % mod;if((k + 2 - i) & 1) ans = (1ll * ans - s[i] + mod) % mod;else ans = (1ll * ans + s[i]) % mod;}printf("%d\n", ans);
}

拉格朗日插值法在ACM竞赛中的一些应用

https://blog.csdn.net/qq_41157137/article/details/100593260?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-1&spm=1001.2101.3001.4242