拉格朗日插值法

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 }

 

 

 

数论的东西，不懂没关系，看还是要看看的...

 

转载于:https://www.cnblogs.com/Judge/p/10428378.html