建议看: cmd:从拉插到快速插值求值
拉格朗日插值
首先我们知道一个事是: $ n+1 $ 个不同点可以唯一确定一个 $ n $ 次多项式 。
那么当你知道了这 $ n+1 $ 个点之后就可以通过拉插插出这个多项式是啥。
定义:
对于 $ n $ 个点 $ (x_i,y_i) $ ,我们希望找到 $ f_i $ 满足:
且 $ f_i $ 也是 $ n-1 $ 次的,那么我们可以得到多项式 $ F = \sum_i f_i $ ,这样 $ F $ 满足他是一个经过这 $ n $ 个点的 $ n-1 $ 次多项式。
那么对于每个 $ f_i $ ,我们尝试用因式构造,简单来说,为了满足 $ f_i(x_j) = 0 ,j \not = i $ ,我们使他含有 $ (x-x_j) , j \not = i $ 的因式,此时再去满足 $ f_i(x_i) = y_i $ ,我们先把 $ x_i $ 带进去,就会得到 $ \prod_{j,j \not = i} (x_i - x_j) $ ,这时我们给整个式子乘上 $ \large \frac{y_i}{\prod_{j,j \not = i} (x_i - x_j)} $ ,就满足了所有条件。
所以我们获得了拉插公式:
例题:
P4781 【模板】拉格朗日插值
码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
typedef unsigned long long ull;
#define mk make_pair
#define ps push_back
#define fi first
#define se second
const int N=1e6+10,inf=0x3f3f3f3f;
const ll mod=998244353,linf=0x3f3f3f3f3f3f3f3f;
inline ll read(){char c=getchar();ll x=0;bool f=0;while(!isdigit(c))f=c=='-'?1:0,c=getchar();while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();return f?-x:x;
}
inline ll qpow(ll x,ll y){ll ans=1;while(y){if(y&1)ans=ans*x%mod;x=x*x%mod;y>>=1;}return ans;
}
int n,K,x[N],y[N];
signed main(){#ifndef ONLINE_JUDGEfreopen("in.in","r",stdin);freopen("out.out","w",stdout);#endifios::sync_with_stdio(0),cin.tie(0),cout.tie(0);n=read(),K=read();for(int i=1;i<=n;++i)x[i]=read(),y[i]=read();ll ans=0;for(int i=1;i<=n;++i){ll i1=1,i2=1;for(int j=1;j<=n;++j){if(j==i)continue;i1=i1*(K-x[j])%mod;i2=i2*(x[i]-x[j])%mod;}i1=(i1+mod)%mod;i2=qpow((i2+mod)%mod,mod-2);ans+=i1*i2%mod*y[i]%mod;}cout<<(ans%mod+mod)%mod;
}
$ x_i $ 取值连续时的插值
我们观察式子
对于 $ \prod x - x_j $ ,我们做一个前缀积和后缀积就可以 $ O(n) $ 解决,对于 $ \prod x_i - x_j $ 因为他们取值连续所以就变成了 $ k_1 ! k_2 ! $ 的形式,同样可以 $ O(n) $ 做。
例题:
CF622F The Sum of the k-th Powers
码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define int ll
typedef pair<int,int> pii;
typedef unsigned long long ull;
#define mk make_pair
#define ps push_back
#define fi first
#define se second
const int N=1e6+10,inf=0x3f3f3f3f;
const ll mod=1e9+7,linf=0x3f3f3f3f3f3f3f3f;
inline ll read(){char c=getchar();ll x=0;bool f=0;while(!isdigit(c))f=c=='-'?1:0,c=getchar();while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();return f?-x:x;
}
inline ll qpow(ll x,ll y){ll ans=1;while(y){if(y&1)ans=ans*x%mod;x=x*x%mod;y>>=1;}return ans;
}
int y[N],n,K,pre[N],hou[N],jc[N];
signed main(){#ifndef ONLINE_JUDGEfreopen("in.in","r",stdin);freopen("out.out","w",stdout);#endifios::sync_with_stdio(0),cin.tie(0),cout.tie(0);n=read(),K=read();y[0]=0;jc[0]=1;for(int i=1;i<=K+2;++i){y[i]=(y[i-1]+qpow(i,K))%mod;jc[i]=jc[i-1]*i%mod;}if(n<=K+2)return cout<<y[n],0;pre[0]=n;for(int i=1;i<=K+2;++i){pre[i]=pre[i-1]*(n-i)%mod;}hou[K+3]=1;for(int i=K+2;i>=0;--i)hou[i]=hou[i+1]*(n-i)%mod;int ans=0;for(int i=0;i<=K+2;++i){ans+=(i==0?1ll:pre[i-1])*hou[i+1]%mod*qpow(jc[i],mod-2)%mod*qpow(jc[K+2-i],mod-2)%mod*y[i]%mod*((K+2-i)&1?-1:1);}cout<<(ans%mod+mod)%mod;
}
重心拉格朗日插值法
感觉没什么意义,就是瞪了两眼式子省去了不必要的计算。
还是看式子:
我们每加入一个点对于每个 $ i $ , $ \prod x - x_j $ 变化可以 $ O(1) $ 算, $ \prod x_i - x_j $ 也可以 $ O(1) $ 算。
插多项式系数
直接插一个值出来是 $ O(n^2) $ 的,但是拉格朗日插值同样可以在 $ O(n^2) $ 的时间复杂度内插出这个多项式的系数,而且支持取模。
还是得改写式子:
首先对于整个式子
有 $ y_i \prod_{j \not = i} {\frac{1}{x_i - x_j} } $ 是常数,姑且写为 $ c_i $ ,可以 $ O(n^2) $ 求出。对于 $ \prod_{j \not = i} x - x_j $ ,写为 $ \frac { \prod_j x - x_j } { x - x_i } $ ,我们把 $ \prod_j x - x_j $ 写为 $ m $ 。
现在式子长成了
首先 $ c_i $ 是常数,不会对 $ x $ 次数造成影响,主要造成影响的是 $ m $ ,它包含 $ n $ 个形如 $ x - x_j $ 的式子,最后又除了一个 $ x - x_i $ ,学过二项式定理的应该知道 $ x $ 的次数就是由每个 $ x - x_j $ 的式子提供来的,考虑 dp ,设 $ dp_i $ 表示 $ x_i $ 的系数,那么每次考虑一个 $ x - x_j $ 的时候,要么选择 $ x $ ,有 $ dp_i \gets dp_{i-1} $ ,如果选择 $ x - x_i $ 的话,就有 $ dp_i = - x_i dp_i $ ,所以最终有 $ dp_i = dp_{i-1} - x_i dp_i $ 。
如果只考虑 $ m $ 的话,直接 dp 是 $ O(n^2) $ 的,但是有了 $ \frac{1}{x - x_i} $ 之后朴素实现就是 $ O(n^3) $ 的了,考虑有 $ \frac{1}{x - x_i} $ 之后相当于一个退背包的过程,所以可以对于每个 $ i $ 直接 $ O(n) $ 做,所以也是 $ O(n^2) $ 的。
码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
#define int ll
typedef pair<int,int> pii;
typedef unsigned long long ull;
#define mk make_pair
#define ps push_back
#define fi first
#define se second
const int N=1e6+10,inf=0x3f3f3f3f;
const ll mod=1e9+7,linf=0x3f3f3f3f3f3f3f3f;
inline ll read(){char c=getchar();ll x=0;bool f=0;while(!isdigit(c))f=c=='-'?1:0,c=getchar();while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar();return f?-x:x;
}
inline ll qpow(ll x,ll y){ll ans=1;while(y){if(y&1)ans=ans*x%mod;x=x*x%mod;y>>=1;}return ans;
}
int n,x[N],y[N],c[N],f[N],g[N],ff[N];
signed main(){#ifndef ONLINE_JUDGEfreopen("in.in","r",stdin);freopen("out.out","w",stdout);#endifios::sync_with_stdio(0),cin.tie(0),cout.tie(0);n=read();for(int i=1;i<=n;++i)x[i]=read(),y[i]=read();for(int i=1;i<=n;++i){c[i]=1;for(int j=1;j<=n;++j)if(j!=i)c[i]=c[i]*(x[i]-x[j])%mod;c[i]=(qpow(c[i],mod-2)*y[i]%mod+mod)%mod;}f[0]=1;for(int i=1;i<=n;++i){for(int j=n-1;j;--j)f[j]=((f[j-1]+f[j]*(-x[i]))%mod+mod)%mod;f[0]=(f[0]*(-x[i])%mod+mod)%mod;}for(int i=1;i<=n;++i){ll ny=qpow(((-x[i])%mod+mod)%mod,mod-2);g[0]=f[0]*ny%mod;ff[0]=(ff[0]+c[i]*g[0])%mod;for(int j=1;j<n;++j){g[j]=((f[j]-g[j-1])*ny%mod+mod)%mod;ff[j]=(ff[j]+c[i]*g[j])%mod;}}for(int i=0;i<n;++i)cout<<(ff[i]+mod)%mod<<' ';
}