拉格朗日插值 学习笔记

概述

\(n\) 个 \(x\) 坐标不同的点可以唯一确定一个 \(n-1\) 次多项式。现在给出了这些点，要求多项式在 \(x\) 处的取值。

考虑构造这个多项式：

首先构造出一个式子，在 \(x_i\) 处有值，而在其他的 \(x_j(j\ne i)\) 处都是零。也就是 \(\prod_{j\ne i}(x-x_j)\)。

此时 \(x_i\) 处的取值是 \(\prod_{j\ne i}(x_i-x_j)\)。乘上一个系数让 \(x_i\) 处的取值为 \(y_i\)，也就是 \(y_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\)。

最后将 \(n\) 个这样的式子加起来就得到了拉插的式子：\(\sum_{i=1}^ny_i\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\)，然后代入 \(x\) 即可。

for(int i=1;i<=n;i++){
  int u=y[i],d=1;
  for(int j=1;j<=n;j++)if(j!=i)u=1ll*u*(k-x[j]+mod)%mod,d=1ll*d*(x[i]-x[j]+mod)%mod;
  ans=(ans+1ll*u*inv(d,mod))%mod;
}

线性插值

有时我们可以选定这 \(n\) 个点的位置，可以选择横坐标 \(1\sim n\)，就可以 \(O(n)\) 插值。

\[\begin{aligned} \sum_{i=1}^ny_i\prod_{j\ne i}\frac{x-j}{x_i-j}&=\sum_{i=1}^ny_i\frac{\prod_{j=1}^{i-1}(x-j)\prod_{j=i+1}^n(x-j)}{\prod_{j=1}^{i-1}(i-j)\prod_{j=i+1}^n(i-j)}\\ &=\sum_{i=1}^ny_i\frac{\prod_{j=1}^{i-1}(x-j)\prod_{j=i+1}^n(x-j)}{(i-1)!(n-i)!(-1)^{n-i}} \end{aligned}\]

预处理 \(x-j\) 的前后缀和与阶乘逆元即可。

CF622F

自然数幂和问题。观察到答案是一个 \(k+1\) 次多项式，可以用上面的方式插 \(k+2\) 个点。

#include<bits/stdc++.h>
using namespace std;
const int mod=1000000007;
int n,k,prime[1000005],y[1000005],fac[1000005],vac[1000005],pre[1000005],suf[1000005],ans;
bool a[1000005];
template<typename T>T qpow(T a,T b,T n,T ans=1){
  for(a%=n;b;b>>=1)b&1&&(ans=1ll*ans*a%n),a=1ll*a*a%n;
  return ans;
}
template<typename T>T inv(T a,T b){
  return qpow(a,b-2,b);
}
int f_init(int n,int prime[],int f[],bool a[],int cnt=0){
  f[1]=1; 
  for(int i=2;i<=n;i++)a[i]=1;
  for(int i=2;i<=n;i++){
    if(a[i])prime[++cnt]=i,f[i]=qpow(i,k,mod);
    for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
      a[i*prime[j]]=0,f[i*prime[j]]=1ll*f[i]*f[prime[j]]%mod;
      if(i%prime[j]==0)break;
    }
  }
  return cnt;
}
int main(){
  cin>>n>>k,f_init(k+2,prime,y,a);
  for(int i=1;i<=k+2;i++)y[i]=(y[i-1]+y[i])%mod; 
  fac[0]=pre[0]=suf[k+3]=1;
  for(int i=1;i<=k+2;i++)pre[i]=1ll*pre[i-1]*(n-i+mod)%mod;
  for(int i=k+2;i>=1;i--)suf[i]=1ll*suf[i+1]*(n-i+mod)%mod;
  for(int i=1;i<=k+2;i++)fac[i]=1ll*fac[i-1]*i%mod;
  vac[0]=1,vac[k+2]=inv(fac[k+2],mod);
  for(int i=k+1;i>=1;i--)vac[i]=1ll*vac[i+1]*(i+1)%mod;
  for(int i=1;i<=k+2;i++)ans=(ans+(((k+2-i)&1?-1ll:1ll)*y[i]*pre[i-1]%mod*suf[i+1]%mod*vac[i-1]%mod*vac[k+2-i]%mod+mod)%mod)%mod;
  return cout<<ans<<'\n',0;
}

[[多项式]]