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 }
数论的东西,不懂没关系,看还是要看看的...