拉格朗日插值法

https://www.cnblogs.com/zwfymqz/p/10063039.html

 

以爲把zwfymqz大佬的博客粘上來就差很少了

 

本博客比較淺顯，適合入門粗學，具體深刻的話就看 attack 大佬的博客（就是上面的連接）吧

 

拉格朗日的公式

 

首先拉格朗日插值法的公式附上：

 

$$A(k)=\sum_{i=1}^{n} y_i \prod_{j=1}^{n} \dfrac{k-x_j}{x_i-x_j}$$

 

這個式子給出來確定是先懵。不要緊，誰看到這玩意兒不會懵呢？

 

拉格朗日的做用

 

解釋一下，拉格朗日插值就是給你 n+1 個點，而後讓你構造出一個符合這堆點練成的光滑圖像的 n 次函數

 

說人話，就是給你 n+1 個點，讓你構造出一個 n 次函數，使得這個函數的圖像通過座標軸上的 n+1 個點。

 

那麼上面式子裏面的 $x_i$ 以及 $y_i$ 就是一個點的座標啦！

 

拉格朗日的分析

 

而後你本身驗證一下就能夠發現： 將原來的點帶進去， A 函數輸出的結果是符合這 n 個點的位置的。

 

爲何？ 由於若是你帶的點是 i 的話，那麼當 $\sum$ 作到第 i 次時，後面的 $\prod$ 算出來的是 1 ，$\sum$的其他項都是 0

 

因而原式成立了，咱們也能夠拿它來計算別的 k 的函數值了

 

 

拉格朗日題的套路

 

至於這玩意兒拿來出題，無非就是裸題：

 

　　給你 n 個點以及正整數 k ，讓你用 n 個點求一個函數圖像並將 k 帶入求值（固然不是真要你求出函數圖像，輸出答案就行了）

 

這種題目的作法要麼就是暴力代入公式，$O(n^2)$ ，要麼給你的 n 個點的橫座標是連續的，你能夠利用這個性質優化到 $O(n)$ ，至於怎麼優化嘛，想一想階乘之類的東西...（具體推導你能夠考慮打開文頭的連接）

 

至於另外一種題型就是讓你求 $\sum_{i=1}^{n}i^k$ 的值了，這個東西爲何能用拉格朗日？他和上面講的東西毫無關聯啊！（你問我我問誰...）總之記好這玩意兒能夠表示成一個 k+1 次多項式，而後就用 n 的值拿進去代就行了

 

作法麼，天然就是枚舉出前 k+2 個點，而後看看 n 是否大於 k+2，不大於的話直接出解了，大於的話就是用 k+2 個點跑拉格朗日，n 帶入就行了。 而且因爲這 k+2 個點的橫座標是連續的，你能夠將計算優化到 $O(k)$

 

什麼？爲何這 k 個點是連續的？ woc 這幾個點是你本身選出來的前 k+2 個點啊！

 

但最後仍是決定解釋一下爲何 $A(n)=\sum_{i=1}^{n}i^k$ 這玩意兒能夠表示成一個 k+1 次函數

 

怎麼說，咱們考慮將 $i^k$ 表示爲 $n-(n-i)^k$ ，懂了吧，將 $n-i$ 看作常數（何況 $n-i$ 在這個表達式中原本就是連續的數）而後展開一下就是一個 k 次多項式啦！

 

恩？ k 次？ 不是 k+1 次麼？ emmm 咱們考慮一下這裏 $n^k$  的項其實總共有 n 項，加起來可不就是 $n*n^k=n^{k+1}$ 麼？ 至於剩下的項次也是能夠以此類推的...

 

（有興趣的讀者能夠算一下哦，畢竟博主沒有算過，算出來piapia 打臉的話還請讀者在評論區中指出，但就算剩下的項沒法整合成高一階的項，這裏多項式的最高次係數也已是 k+1 ，知足 k+1 次函數的定義了）

 

神奇的重心拉格朗日

 

而後是某種騷操做，我是搞不大懂，拉格朗日前面加了個重心，說是優化，而後又說每加一個點的複雜度是 $O (n)$ 的（那加 n 個點的複雜度不仍是 $n^2$ ?），也不知道真的假的...貌似真的

 

具體到圖像上的實現的話有一個動圖能夠說明（不知道放到cnblogs上能不能動起來）：

 

 

 

 

怎麼說呢...就是慢慢找點唄...

 

而後就 完 了 _(:з」∠)_

 

慢着？板子都不放一個的？

 

$$O(n^2)$$

 1 //by Judge
 2 //by young Judge
 3 #include<iostream>
 4 #include<cstdio>
 5 #define ll long long
 6 using namespace std;
 7 const int M=2111;
 8 const ll mod=998244353;
 9 #ifndef Judge
10 #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
11 #endif
12 char buf[1<<21],*p1=buf,*p2=buf;
13 inline ll read(){
14     ll x=0,f=1; char c=getchar();
15     for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
16     for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
17 } int n; ll k,K,ans,x[M],y[M];
18 inline void AMOD(ll& a,ll& b){ a+=b; if(a>=mod) a%=mod; }
19 inline ll quick_pow(ll x,ll p,ll ans=1){
20     while(p){
21         if(p&1) ans=ans*x%mod;
22         x=x*x%mod,p>>=1;
23     } return ans;
24 }
25 int main(){
26     n=read(),K=read();
27     for(int i=1;i<=n;++i) x[i]=read(),y[i]=read();
28     for(int i=1,j;i<=n;++i){ k=1;
29         for(j=1;j<=n;++j) if(i!=j)
30             k=k*(x[i]+mod-x[j])%mod;
31         k=quick_pow(k,mod-2);
32         for(j=1;j<=n;++j) if(i!=j)
33             k=k*(K+mod-x[j])%mod;
34         k=k*y[i]%mod,AMOD(ans,k);
35     } return printf("%lld\n",ans),0;
36 }

 

$$O(n)$$

 1 //by Judge
 2 #include<iostream>
 3 #include<cstdio>
 4 #define ll long long
 5 using namespace std;
 6 const int M=2111;
 7 const ll mod=998244353;
 8 #ifndef Judge
 9 #define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
10 #endif
11 char buf[1<<21],*p1=buf,*p2=buf;
12 inline ll read(){
13     ll x=0,f=1; char c=getchar();
14     for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
15     for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
16 } ll n,K,ans,x[M],y[M],s1[M],s2[M],ifac[M];
17 inline void ADD(ll& a,ll b){ a+=b; if(a>=mod) a%=mod; }
18 int main(){ n=read()-1,K=read(),s2[n+1]=ifac[0]=ifac[1]=1;
19     for(int i=0;i<=n;++i) x[i]=read(),y[i]=read(); s1[0]=(K-x[i])%mod;
20     for(int i=1;i<=n;++i) s1[i]=s1[i-1]*(K-x[i])%mod;
21     for(int i=n;i>=1;--i) s2[i]=s2[i+1]*(K-x[i])%mod;
22     for(int i=2;i<=n;++i) ifac[i]=(mod-mod/i)*ifac[mod%i]%mod;
23     for(int i=2;i<=n;++i) ifac[i]=ifac[i]*ifac[i-1]%mod;
24     for(int i=0;i<=n;++i)
25         ADD(ans,1ll*y[i]*(i?s1[i-1]:1)%mod*s2[i+1]%mod*
26             ifac[i]%mod*((n-i&1)?-1:1)*ifac[n-i]%mod);
27     return !printf("%lld\n",(ans+mod)%mod);
28 }

 

什麼？要函數封裝的？

 1 ll lagrange(ll n,ll *x,ll *y,ll xi){ ll ans=0;
 2     s1[0]=(xi-x[0])mod,s2[n+1]=ifac[0]=ifac[1]=1;
 3     for(int i=1;i<=n;++i) s1[i]=s1[i-1]*(xi-x[i])%mod;
 4     for(int i=n;i>=0;--i) s2[i]=s2[i+1]*(xi-x[i])%mod;
 5     for(int i=2;i<=n;++i) ifac[i]=(mod-mod/i)*ifac[mod%i]%mod;
 6     for(int i=2;i<=n;++i) ifac[i]=ifac[i]*ifac[i-1]%mod;
 7     for(int i=0;i<=n;++i)
 8         ADD(ans,1ll*y[i]*(i?s1[i-1]:1)%mod*s2[i+1]%mod*
 9             ifac[i]%mod*((n-i&1)?-1:1)*ifac[n-i]%mod);
10     return (ans+mod)%mod;
11 }

 

 

 

數論的東西，不懂不要緊，看仍是要看看的...