zoukankan      html  css  js  c++  java
  • 拉格朗日插值法

    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上能不能动起来):

    重心拉格朗日gif

    怎么说呢...就是慢慢找点呗...

    然后就 完 了 _(:з」∠)_

    慢着?板子都不放一个的?

    $$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
    ",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
    ",(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 }

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

  • 相关阅读:
    Linux内核初探 之 进程(三) —— 进程调度算法
    Android中的路径记录 | RobinBlog
    子域名劫持
    zookeeper 实战
    [iOS 开发] WebViewJavascriptBridge 从原理到实战 · Shannon's Blog
    Swift Property
    工厂方法模式
    jquery插件封装
    其他事件
    吴裕雄--天生自然诗经学习笔记 :夸父逐日
  • 原文地址:https://www.cnblogs.com/Judge/p/10428378.html
Copyright © 2011-2022 走看看