- C20250053's blog
拉格朗日插值学习笔记
- @ 2026-10-5 15:06:56
拉格朗日插值
给出 oi-wiki 中插值的定义:
插值是一种通过已知的、离散的数据点推算一定范围内的新数据点的方法。插值法常用于函数拟合中。
拉格朗日插值是一种多项式插值。
给定 个点 ,若满足 ,那么这 个点可以唯一确定一个 次多项式 ,求这个多项式 ,并给定参数 ,求 的值。
我们考虑 构造 满足:
$$f_i(x_j)=\begin{cases}y_j,j=i\\ 0,j\not=i \end{cases}$$此时 ,满足 。
我们尝试构造 :
$$f_i(x)=\frac{y_i}{\prod_{j\not=i}(x_i-x_j)}\times \prod_{j\not=i}(x-x_j)\\ f_i(x)=y_i\times \frac{\prod_{j\not=i}(x-x_j)}{\prod_{j\not=i}(x_i-x_j)}\\ =y_i\times \prod_{j\not=i}\frac{x-x_j}{x_i-x_j}$$所以 $f(x)=\displaystyle\sum_{i=1}^n y_i\times \prod_{j\not=i}\dfrac{x-x_j}{x_i-x_j}$。
附上模板题代码:
#include<bits/stdc++.h>
#define FL(i,a,b) for(int i=(a);i<=(b);i++)
#define FR(i,a,b) for(int i=(a);i>=(b);i--)
#define ll long long
#define ull unsigned long long
#define ld long double
#define PII pair<int,int>
using namespace std;
const int MAXN = 2e3 + 10;
const int mod = 998244353;
int n,K,ans=0;
int x[MAXN],y[MAXN];
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 main(){
scanf("%d%d",&n,&K);
FL(i,1,n) scanf("%d%d",&x[i],&y[i]);
FL(i,1,n){
ll p=y[i],q=1;
FL(j,1,n){
if(i!=j){
p=1ll*p*(K-x[j])%mod;
q=1ll*q*(x[i]-x[j])%mod;
}
}
ans=(ans+1ll*p*qpow(q,mod-2)%mod)%mod;
}
printf("%d\n",(ans+mod)%mod);
}
CF1874E Jellyfish and Hack
感觉这道题的 trick 还是很新奇的。
考虑 dp,设 表示长度为 且 的排列 的数量。
枚举 ,把序列划分为长度 和 的两个子序列(即 是 序列中第 大的数)。
发现合法的 的量级是 的。
转移方程如下:
$dp_{x,y}=\sum_{i=1}^x\binom{x-1}{i-1}\sum_{j=0}^{y-x}dp_{i-1,j}\times dp_{x-i,y-x-j}$。
朴素转移是 的,考虑优化。
我们令 $F_n(t)=\displaystyle \sum_{i=0}^{\infty}dp_{n,i}t^i$,则有:
$$F_n(t)=\sum_{v=0}^{\infty}\sum_{i=1}^n\binom{n-1}{i-1}\sum_{j=0}^{v-n}dp_{i-1,j}\times dp_{n-i,v-n-j}\times t^v\\ F_n(t)=\sum_{i=1}^n\binom{n-1}{i-1}\sum_{v=0}^{\infty}\sum_{j=0}^{v-n}dp_{i-1,j}\times dp_{n-i,v-n-j}\times t^v\\ 令 k=v-n,\\ F_n(t)=\sum_{i=1}^n\binom{n-1}{i-1}\sum_{v=0}^{\infty}\sum_{j=0}^{k}dp_{i-1,j}\times dp_{n-i,k-j}\times t^{n+k}\\ F_n(t)=\sum_{i=1}^n\binom{n-1}{i-1}\sum_{k=0}^{\infty}(\sum_{j=0}^{k}dp_{i-1,j}\times dp_{n-i,k-j}\times t^k)t^n\\ F_n(t)=\sum_{i=1}^n\binom{n-1}{i-1}F_{i-1}(t)F_{n-i}(t)t^n$$以及多项式乘法本质上就是序列的卷积。
然后这个多项式的 是 的(多项式的度()是指该多项式中最高次项的次数)
我们考虑带入 个值来求出这个多项式的系数。
我们枚举 ,预处理 的幂次,可以 算出 。
我们回忆拉格朗日插值的公式,然后对于每个 代入 ,用拉格朗日插值可以 算出 的各项系数,即 。
$F_n(x)=\displaystyle\sum_{i=1}^{n^2+\varepsilon} y_i\times \prod_{j\not=i}\dfrac{x-x_j}{x_i-x_j}$。
我们令 $a_i=y_i\times \displaystyle\prod_{j\not=i} \frac{1}{x_i-x_j},L(x)=\prod_{j=1}^n(x-x_j)$,则 $F_n(x)=\displaystyle\sum_{i=1}^{n^2+\varepsilon} a_i\times \frac{L(x)}{x-x_i}$。
我们计算 中 次项的系数 ,然后对每个数 ,做多项式除法(除以 ,然后因为一定是能整除的所以可以正着做)求出 次项的系数 ,于是 第 项的系数(即 )加上 。
最后答案即为 $\displaystyle \sum_{i=lim}^{n^2+\varepsilon}dp_{n,i}$。
#include<bits/stdc++.h>
#define FL(i,a,b) for(int i=(a);i<=(b);i++)
#define FR(i,a,b) for(int i=(a);i>=(b);i--)
#define ll long long
#define ull unsigned long long
#define ld long double
#define PII pair<int,int>
using namespace std;
const int MAXN = 2e2 + 10;
const int mod = 1e9 + 7;
const int eplison = 5;
int n,lim,Mx,ans=0;
int C[MAXN][MAXN],pw[MAXN*MAXN][MAXN];
int x[MAXN*MAXN],y[MAXN*MAXN];
int F[MAXN][MAXN*MAXN],cnt[MAXN*MAXN];
int a[MAXN*MAXN],s[MAXN*MAXN],b[MAXN*MAXN];
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;
}
void Solve(){
FL(i,1,Mx) x[i]=i,y[i]=F[n][i];
FL(i,1,Mx){
int Mul=1;
FL(j,1,Mx)
if(i!=j) Mul=1ll*Mul*(x[i]-x[j]+mod)%mod;
a[i]=1ll*y[i]*qpow(Mul,mod-2)%mod;
}
s[0]=1;
FL(i,1,Mx){
FR(j,i,1)
s[j]=(s[j-1]-1ll*x[i]*s[j]%mod+mod)%mod;
s[0]=1ll*(mod-s[0])*x[i]%mod;
}
FL(i,1,Mx){
int inv=qpow(mod-x[i],mod-2);
b[0]=1ll*s[0]*inv%mod;
cnt[0]=(cnt[0]+1ll*a[i]*b[0]%mod)%mod;
FL(j,1,Mx){
b[j]=1ll*(s[j]-b[j-1]+mod)*inv%mod;
cnt[j]=(cnt[j]+1ll*a[i]*b[j]%mod)%mod;
}
}
}
int main(){
scanf("%d%d",&n,&lim),Mx=n*(n+1)/2;
if(lim>Mx) puts("0"),exit(0);
Mx+=eplison;
FL(i,1,Mx){
pw[i][0]=1;
FL(j,1,n) pw[i][j]=1ll*pw[i][j-1]*i%mod;
}
FL(i,0,MAXN-1){
C[i][0]=C[i][i]=1;
FL(j,1,i-1) C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
}
FL(i,1,Mx) F[1][i]=i,F[0][i]=1;
FL(i,2,n){
FL(j,1,Mx){
FL(k,1,i)
F[i][j]=(F[i][j]+1ll*C[i-1][k-1]*F[k-1][j]%mod*F[i-k][j]%mod)%mod;
F[i][j]=1ll*F[i][j]*pw[j][i]%mod;
}
}
Solve();
FL(i,lim,Mx) ans=(ans+cnt[i])%mod;
printf("%d\n",ans);
}
省集题目(无题号)
给定 个节点的树,每个节点有权值 ,终点在树的一个叶子节点,但是对你未知。
你需要在走到终点的时候恰好选择 个节点,且选择的节点的 按选择顺序单调不增。
但是对于每个节点 ,只有在你对其进行是否选择的决策后你才能得知终点在 的哪个儿子的子树中。
所以你会选择一个节点集合 ,在经过集合中的节点时必然选择这个节点,在经过集合外的节点时必然不选择这个节点。
并且保证每个叶子节点作为终点的情况都能满足条件。
问对于所有可能的权值序列 ,选取节点集合的方案数的总和,对 取模。
。
我们称在集合 中的节点为关键点,要求根到每个叶子的路径都满足恰好经过 个关键点,且关键点权值按路径顺序单调不增。
我们考虑朴素树形 dp,令 表示在 子树内给每个节点赋值并选点,满足 到每个叶子的路径上恰好选 个点,且所有关键点权值 的方案数。
我们令 ,考虑 是否选入集合:
- 如果不选 ,那么 不影响 的值,可以随便填,贡献为 。
- 否则令 ,贡献为 。
因此 的转移式子是:
$$f_{u,i,j}=(R_u-L_u+1)\cdot g_{u,i,j}+\sum_{y=L_u}^{\min(R_u,j)}g_{u,i-1,y}$$时间复杂度 ,可以前缀和优化到 ,其中 。
发现 的取值范围是 的,直接枚举显然不可行,所以考虑离散化,把 放进离散化数组里得到 个不同的数,中间有 段。
发现同一段中的数,和各个节点的取值范围的包含关系不变。
考虑对每段分别处理,设当前在处理第 段,也就是 ,发现 是关于 是一个次数不超过 的多项式 (这个可以从叶子开始推,每个点的选择方案是带 的式子,显然乘起来是多项式),于是我们求出其在 的点值然后用拉格朗日插值还原多项式并求出整段的答案。
具体地,我们固定 ,维护 ,定义 表示在处理当前段 之前,所有已经处理过且属于 的权值范围的 的 的和。
对于节点 :
-
如果 的取值范围覆盖这一段,那么对于上界 ,选 的贡献就是 。
-
否则选 的贡献只有以前贡献和 。
转移后,我们需要把这段的新值加入 :
我们已经得到 时 的值,用拉格朗日插值可以求出 ,然后带入 ,此时多项式等于 $S_{u,i}(p)=his_{u,i}+\sum_{q=0}^{b_{X+1}-b_X-1}g_{u,i,b_X+q}$,此时对应的 ,正好是我们想求的 的新值。
然后最后用同样的方法还原答案,即令 ,求 的值,此时对应的 ,也就是还没有选择任何点且没有任何取值限制时的答案。
时间复杂度 。
#include<bits/stdc++.h>
#define FL(i,a,b) for(int i=(a);i<=(b);i++)
#define FR(i,a,b) for(int i=(a);i>=(b);i--)
#define ll long long
#define ull unsigned long long
#define ld long double
#define PII pair<int,int>
using namespace std;
const int MAXN = 2e2 + 10;
const int MAXK = 20 + 5;
const int mod = 1004535809;
int n,K,X;
int L[MAXN],R[MAXN],sum[MAXN];
int f[MAXN][MAXK][MAXN],g[MAXK][MAXN];
int his[MAXN][MAXK];
vector<int>G[MAXN];
int b[MAXN<<1],tot=0;
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 suf[MAXN],pre[MAXN],exc[MAXN];
void init(){
FL(i,0,n+1){
int Mul=1;
FL(j,0,n+1)
if(j!=i) Mul=1ll*Mul*(i-j+mod)%mod;
exc[i]=qpow(Mul,mod-2);
}
}
int Get(int *f,int ps){
int R=b[ps+1]-b[ps]-1;
suf[n+1]=(R-(n+1)+mod)%mod;
FR(i,n,0) suf[i]=1ll*suf[i+1]*(R-i+mod)%mod;
pre[0]=R;
FL(i,1,n+1) pre[i]=1ll*pre[i-1]*(R-i+mod)%mod;
int Mul=1,res=0;
FL(i,0,n+1){
res=(res+1ll*f[i]*exc[i]%mod*(i?pre[i-1]:1)%mod*(i+1<=n+1?suf[i+1]:1)%mod)%mod;
Mul=1ll*Mul*(R-i+mod)%mod;
}
return res;
}
void dfs1(int u,int fth){
for(int v:G[u]){
if(v==fth) continue;
dfs1(v,u),sum[u]++;
}
}
void dfs2(int u,int fth){
for(int v:G[u]){
if(v==fth) continue;
dfs2(v,u);
}
FL(j,0,K) FL(k,0,n+1) g[j][k]=(!sum[u]&&j?0:1);
for(int v:G[u]){
if(v==fth) continue;
FL(j,0,K) FL(k,0,n+1) g[j][k]=1ll*g[j][k]*f[v][j][k]%mod;
}
int len=b[R[u]+1]-b[L[u]];
FL(j,0,K) FL(k,0,n+1) f[u][j][k]=1ll*len*g[j][k]%mod;
FL(j,0,K-1){
int prv=his[u][j];
if(L[u]<=X&&X<=R[u]){
int sum=prv;
FL(k,0,n+1){
sum=(sum+g[j][k])%mod;
f[u][j+1][k]=(f[u][j+1][k]+sum)%mod,g[j][k]=sum;
}
his[u][j]=Get(g[j],X);
}
else FL(k,0,n+1) f[u][j+1][k]=(f[u][j+1][k]+prv)%mod;
}
}
int main(){
scanf("%d%d",&n,&K);
FL(i,1,n) scanf("%d",&L[i]);
FL(i,1,n) scanf("%d",&R[i]);
FL(i,1,n-1){
int u,v;
scanf("%d%d",&u,&v);
G[u].push_back(v);
G[v].push_back(u);
}
FL(i,1,n) b[++tot]=L[i],b[++tot]=R[i]+1;
sort(b+1,b+tot+1);
tot=unique(b+1,b+tot+1)-b-1;
FL(i,1,n)
L[i]=lower_bound(b+1,b+tot+1,L[i])-b,
R[i]=lower_bound(b+1,b+tot+1,R[i]+1)-b-1;
init();
dfs1(1,0);
FL(i,1,tot-1) X=i,dfs2(1,0);
printf("%d\n",Get(f[1][K],tot-1));
}