文章目錄

• • 介紹
  • 拉格朗日插值法
  • 重心拉格朗日插值法
  • • 例題1
    • • 技巧
    • 例題2
    • 例題3
  • 練習(xí)題

介紹

眾所周知，$n$ 個點 $(x_i,y_i)$ 可以確定唯一一個 $n?1$ 次的多項式。

拉格朗日插值法就是，用這 $n$ 個點表示出這個多項式。

拉格朗日插值法

考慮構(gòu)造 $n$ 個多項式，第 $i$ 個多項式 $f_i(x)$ 滿足：當 $x=x_i$ 時，$f_i(x)=y_i$，當 $x=x_j(j=i)$ 時，$f_i(x)=0$。

然后將這 $n$ 個多項式加起來，就能得到一個多項式 $f(x)$，滿足當 $x=x_i$ 時，$f(x)=y_i$；并且這是個 $n?1$ 次多項式。

$f_i(x)$ 的構(gòu)造為：$f_i(x)=y_i{\prod }_{j=i}\frac{x?x_j}{x_i?x_j}$。

加起來，得到 $f(x)={\sum }_{i=1}^n{y}_i{\prod }_{j=i}\frac{x?x_j}{x_i?x_j}$。

構(gòu)造多項式 $n^2$，求解一次 $n^2$。

重心拉格朗日插值法

稍微改改上面那個柿子：
$$f(x)=\sum _{i=1}^n{y}_i\prod _{j=i}\frac{x?x_j}{x_i?x_j}$$

設(shè) $h={\prod }_{i=1}^{n}x?x_i$，帶入得：
$$=h\sum _{i=1}^n\prod _{j=i}\frac{y_i}{(x_i?x_j)(x?x_i)}$$

設(shè) $t_i={\prod }_{j=i}\frac{y_i}{x_i?x_j}$，帶入得：
$$=h\sum _{i=1}^n\frac{t_i}{x?x_i}$$

每次加入一個新點時，計算出它的 $t_i$，并且更新別的點的 $t_i$，時間復(fù)雜度 $O(n)$，加入 $n$ 個點就是 $n^2$。

求解時，$O(n)$ 求 $h$，$O(n)$ 求出那個 $\sum$，時間復(fù)雜度 $O(n)$。

例題1

【模板】拉格朗日插值

寫法是重心拉格朗日插值法，雖然在這題會慢，但只是權(quán)當練習(xí)而已，無所謂了。

#include <cstdio> #include <cstring> #include <algorithm> using namespace std; #define mod 998244353int n,k; int ksm(int x,int y){int re=1;for(;(y&1?re=1ll*re*x%mod:0),y;y>>=1,x=1ll*x*x%mod);return re;} #define inv(x) ksm(x,mod-2) struct Lagrange{struct Point{int x,y,t;}d[2010];int tot=0;Lagrange(){}void insert(int x,int y){for(int i=1;i<=tot;i++)d[i].t=1ll*d[i].t*inv(d[i].x-x+mod)%mod;d[++tot]=(Point){x,y,y};for(int i=1;i<tot;i++)d[tot].t=1ll*d[tot].t*inv(x-d[i].x+mod)%mod;}int calc(int x){int prod=1,re=0;for(int i=1;i<=tot;i++){prod=1ll*prod*(x-d[i].x+mod)%mod;re=(re+1ll*d[i].t*inv(x-d[i].x+mod)%mod)%mod;}return 1ll*re*prod%mod;} }Lag;int main() {scanf("%d %d",&n,&k);for(int i=1,x,y;i<=n;i++){scanf("%d %d",&x,&y);Lag.insert(x,y);}printf("%d",Lag.calc(k)); }

技巧

有時候我們遇到的 $n$ 個點的 $x$ 值是連續(xù)的，下面不妨當成 $x_i=i$。

再看上面的柿子：
$$\begin{array}{rl}f(x)& =\sum _{i=1}^n{y}_i\prod _{j=i}\frac{x?x_j}{x_i?x_j}\\ & =\sum _{i=1}^n{y}_i\prod _{j=i}\frac{x?j}{i?j}\end{array}$$

預(yù)處理兩個數(shù)組：$pr{e}_i={\prod }_{j=1}^i(x?j),su{f}_i={\prod }_{j=i}^n(x?j)$，那么分子部分就是 $pr{e}_{j?1}\times su{f}_{j+1}$。

而分母部分不難發(fā)現(xiàn)是個階乘，即 $(i?1)!(n?i)!(?1{)}^{n?i}$。

代入得到：
$$f(x)=\sum _{i=1}^n{y}_i\frac{pr{e}_{i?1}\times su{f}_{i+1}}{(i?1)!(n?i)!(?1{)}^{n?i}}$$

這樣在一頓預(yù)處理后也可以做到 $O(n)$，但他的價值不止如此。

例題2

拉格朗日插值2

當 $i$ 從 $0$ 開始計數(shù)時，上面那個柿子就變成了 $f(x)={\sum }_{i=0}^n{y}_i\frac{pr{e}_{i?1}\times su{f}_{i+1}}{i!(n?i)!(?1{)}^{n?i}}$。

則 $f(m+i)={\prod }_{j=0}^n(m+i?j){\sum }_{j=0}^n\frac{y_j}{j!(n?j)!(?1{)}^{n?j}}\times \frac{1}{m+i?j}$。

前面的 $\prod$ 易于處理，可以先不管，后面的 $\sum$ 是個卷積的形式，但是 $m$ 高達 $10^8$，直接卷不太行。

但實際上處理方法也很簡單，原來是設(shè) $A_i=\frac{y_i}{i!(n?i)!(?1{)}^{n?i}},B_i=\frac{1}{i}$，注意到只有 $[m?n,m+n]$ 區(qū)間內(nèi)的 $B$ 是有用的，不妨將 $B$ 整體左移 $m?n$ 位，即 $B_i=\frac{1}{i+m?n}$，那么 $A$ 卷 $B$ 得到的 $C$ 中，$C_{i+n}$ 即 $f_{m+i}$。

代碼如下：

#include <cstdio> #include <cstring> #include <cmath> #include <algorithm> using namespace std; #define maxn 550010 #define mod 998244353 #define bin(x) (1<<(x))int n,m,A[maxn],B[maxn]; int ksm(int x,int y){int re=1;for(;(y&1?re=1ll*re*x%mod:0),y;y>>=1,x=1ll*x*x%mod);return re;} int fac[maxn],inv_fac[maxn],inv[maxn]; void FacInit(){fac[0]=inv_fac[0]=1;for(int i=1;i<=n;i++)fac[i]=1ll*fac[i-1]*i%mod;for(int i=1;i<=n;i++)inv_fac[i]=1ll*inv_fac[i-1]*inv[i]%mod; } int w[maxn];void prep(int lg){int N=bin(lg);inv[1]=1;for(int i=2;i<=maxn-10;i++)inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;for(int i=1,wn;i<N;i<<=1){w[i]=1;wn=ksm(3,(mod-1)/(i<<1));for(int j=1;j<i;j++)w[i+j]=1ll*w[i+j-1]*wn%mod;} } int r[maxn],limit; void InitR(int lg){for(int i=1;i<bin(lg);i++)r[i]=(r[i>>1]>>1)|((i&1)<<(lg-1));} int add(int x){return x>=mod?x-mod:x;} int dec(int x){return x<0?x+mod:x;} void ntt(int *f,int lg,int type=0){limit=bin(lg);if(type)reverse(f+1,f+limit);for(int i=1;i<limit;i++)if(i<r[i])swap(f[i],f[r[i]]);for(int mid=1,t;mid<limit;mid<<=1)for(int j=0;j<limit;j+=(mid<<1))for(int i=0;i<mid;i++){t=1ll*f[j+i+mid]*w[mid+i]%mod;f[j+i+mid]=dec(f[j+i]-t);f[j+i]=add(f[j+i]+t);}if(type)for(int i=0;i<limit;i++)f[i]=1ll*f[i]*inv[limit]%mod; } void NTT(int *f,int *g,int ln){int lg=ceil(log2(ln));InitR(lg);ntt(f,lg);ntt(g,lg);for(int i=0;i<bin(lg);i++)f[i]=1ll*f[i]*g[i]%mod;ntt(f,lg,1); }int main() {scanf("%d %d",&n,&m);prep(ceil(log2((n+1)<<1)));FacInit();for(int i=0,x;i<=n;i++){scanf("%d",&x);A[i]=1ll*x*inv_fac[i]%mod*inv_fac[n-i]%mod;if(n-i&1)A[i]=(mod-A[i])%mod;}for(int i=0;i<=(n<<1);i++)B[i]=ksm(m-n+i,mod-2);NTT(A,B,2*n+1);int pd=1;for(int i=0;i<=n;i++)pd=1ll*pd*(m-i)%mod;for(int i=n;i<=(n<<1);i++){A[i]=1ll*A[i]*pd%mod;printf("%d ",A[i]);pd=1ll*pd*ksm(m+i-2*n,mod-2)%mod*(m+i-n+1)%mod;} }

例題3

還有一個相當經(jīng)典的應(yīng)用：自然數(shù)冪求和，即給出 $n,k$，求 ${\sum }_{i=1}^n{i}^k$。

這是一個關(guān)于 $n$ 的 $k+1$ 次多項式，證明可以看洛谷第一篇題解，寫的很好，大致思路是：一個 $n$ 次多項式 $f(x)$，假如有一個數(shù)列 $f(1),f(2),f(3),...$，將其差分后得到 $f(2)?f(1),f(3)?f(2),...$ ，這個差分后的數(shù)列的每一項就變成了一個關(guān)于 $i$ 的 $n?1$ 次多項式。而題中的數(shù)列即為 $f(i)={\sum }_{j=1}^i{j}^k$，做差后 $\Delta f(i)=(i+1{)}^k$，是一個關(guān)于 $i$ 的 $k$ 次多項式，所以原數(shù)列是個關(guān)于 $n$ 的 $k+1$ 次多項式。

考慮拉格朗日插值法，那么需要找到 $k+2$ 個點值，由于這里可以自己選取點，所以為了能套用上面的技巧，選取的點就是 $(1,f(1)),(2,f(2)),...,(k+2,f(k+2))$。

代入公式：
$$F(x)=\sum _{i=1}^{k+2}f(i)\frac{pr{e}_{i?1}\times su{f}_{i+1}}{(i?1)!(k+2?i)!(?1{)}^{k+2?i}}$$

再把 $n$ 代入求值就做完了，時間復(fù)雜度 $O(k\log k)$。

當然也可以做到線性：注意到 $i^k$ 是個完全積性函數(shù)，即 $i^k{j}^k=(ij{)}^k$，所以可以線性篩出來，那么就可以線性求 $f$ 了，然后階乘的逆元也是可以 $O(k)$ 預(yù)處理的，這樣就做到了 $O(k)$，但是比較麻煩所以下面的代碼還是寫的 $O(k\log k)$ 做法。

#include <cstdio> #include <cstring> #include <algorithm> using namespace std; #define maxn 1000010 #define mod 1000000007int n,k,ans=0; int ksm(int x,int y){int re=1;for(;(y&1?re=1ll*re*x%mod:0),y;y>>=1,x=1ll*x*x%mod);return re;} void add(int &x,int y){x=(x+y>=mod?x+y-mod:x+y);} void dec(int &x,int y){x=(x-y<0?x-y+mod:x-y);} int pre[maxn],suf[maxn],fac[maxn],inv_fac[maxn];int main() {scanf("%d %d",&n,&k);if(n<=k+2){for(int i=1;i<=n;i++)add(ans,ksm(i,k));}else{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>=1;i--)suf[i]=1ll*suf[i+1]*(n-i)%mod;fac[0]=1;for(int i=1;i<=k+2;i++)fac[i]=1ll*fac[i-1]*i%mod;inv_fac[k+2]=ksm(fac[k+2],mod-2);inv_fac[0]=1;for(int i=k+1;i>=1;i--)inv_fac[i]=1ll*inv_fac[i+1]*(i+1)%mod;int sum=0;for(int i=1;i<=k+2;i++){add(sum,ksm(i,k));(k+2-i&1?dec:add)(ans,1ll*sum*pre[i-1]%mod*suf[i+1]%mod*inv_fac[i-1]%mod*inv_fac[k+2-i]%mod);}}printf("%d",ans); }

練習(xí)題

都是些眾所周知的題，以后要是遇到什么有趣的也會放這兒的。

[TJOI2018]教科書般的褻瀆 ? 題解
tyvj 1858 XLkxc ? 題解
[NOIP2020] 微信步數(shù) ? 題解
[集訓(xùn)隊互測2012] calc ? 題解
[JLOI2016]成績比較 ? 題解

總結(jié)

以上是生活随笔為你收集整理的拉格朗日插值法小结的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。