拉格朗日插值法及应用

拉格朗日插值法

  • 快速根据点值逼近函数
  • 在取点大于 n n <script type="math/tex" id="MathJax-Element-12">n</script>的情况下解出 n n <script type="math/tex" id="MathJax-Element-13">n</script>次多项式是唯一解。

例如:
i=1ni=n(n+1)2 ∑ i = 1 n i = n ( n + 1 ) 2 <script type="math/tex" id="MathJax-Element-14">\sum\limits_{i=1}^ni=\frac{n(n+1)}{2}</script>
只需取三个点 (1,1),(2,3),(3,6) ( 1 , 1 ) , ( 2 , 3 ) , ( 3 , 6 ) <script type="math/tex" id="MathJax-Element-15">(1,1),(2,3),(3,6)</script>即可确定唯一方程。

同理 i=1ni2 ∑ i = 1 n i 2 <script type="math/tex" id="MathJax-Element-16">\sum\limits_{i=1}^{n}i^2</script>只需由四个点确定,但可以取三个点来逼近。

作用?
当n很大时,直接带入求出的函数求解。

一般方法

对于 n+1 n + 1 <script type="math/tex" id="MathJax-Element-17">n+1</script>个点对 (x0,y0),(x1,y1)...(xn,yn) ( x 0 , y 0 ) , ( x 1 , y 1 ) . . . ( x n , y n ) <script type="math/tex" id="MathJax-Element-18">(x_0,y_0),(x_1,y_1)...(x_n,y_n)</script>,求一个函数 fi f i <script type="math/tex" id="MathJax-Element-19">f_i</script>,使得该函数在 xi x i <script type="math/tex" id="MathJax-Element-20">x_i</script>处取得对应 yi y i <script type="math/tex" id="MathJax-Element-21">y_i</script>值,而在其他 xj x j <script type="math/tex" id="MathJax-Element-22">x_j</script>取得 0 0 <script type="math/tex" id="MathJax-Element-23">0</script>,最后把这几个函数线性结合即可。可得:

fi(x)=ji(xxj)ji(xixj)yi f i ( x ) = ∏ j ≠ i ( x − x j ) ∏ j ≠ i ( x i − x j ) ∗ y i
<script type="math/tex; mode=display" id="MathJax-Element-24">f_i(x)=\frac{\prod\limits_{j\not = i}(x-x_j)}{\prod\limits_{j\not = i}(x_i-x_j)}*y_i</script>

g(x)=i=0nfi(x) g ( x ) = ∑ i = 0 n f i ( x )
<script type="math/tex; mode=display" id="MathJax-Element-25">g(x)=\sum_{i=0}^nf_i(x)</script>

求原函数的复杂度为 O(n2) O ( n 2 ) <script type="math/tex" id="MathJax-Element-26">O(n^2)</script>,单次求值的复杂度为 O(n2) O ( n 2 ) <script type="math/tex" id="MathJax-Element-27">O(n^2)</script>。

重心拉格朗日插值法

考虑得到的函数:

fi(x)=ji(xxj)ji(xixj)yi f i ( x ) = ∏ j ≠ i ( x − x j ) ∏ j ≠ i ( x i − x j ) ∗ y i
<script type="math/tex; mode=display" id="MathJax-Element-28">f_i(x)=\frac{\prod\limits_{j\not = i}(x-x_j)}{\prod\limits_{j\not = i}(x_i-x_j)}*y_i</script>

g(x)=i=0nfi(x) g ( x ) = ∑ i = 0 n f i ( x )
<script type="math/tex; mode=display" id="MathJax-Element-84">g(x)=\sum_{i=0}^nf_i(x)</script>

不难发现每个 fi f i <script type="math/tex" id="MathJax-Element-85">f_i</script>计算了重复部分。

l(x)=i=1n(xxi) l ( x ) = ∏ i = 1 n ( x − x i )
<script type="math/tex; mode=display" id="MathJax-Element-86">l(x) = \prod_{i=1}^{n}(x-x_i)</script>
对于每一个 fi f i <script type="math/tex" id="MathJax-Element-87">f_i</script>
设:
ωi=yiji(xixj) ω i = y i ∏ j ≠ i ( x i − x j )
<script type="math/tex; mode=display" id="MathJax-Element-88">\omega_i= \frac{y_i}{\prod\limits_{j\not = i}(x_i-x_j)}</script>

则:

g(x)=l(x)i=1nωi(xxi) g ( x ) = l ( x ) ∑ i = 1 n ω i ( x − x i )
<script type="math/tex; mode=display" id="MathJax-Element-89">g(x)=l(x)\sum_{i=1}^{n}\frac{\omega _i}{(x-x_i)}</script>

求原函数复杂度为 O(n2) O ( n 2 ) <script type="math/tex" id="MathJax-Element-90">O(n^2)</script>,单次求值复杂度 O(n2) O ( n 2 ) <script type="math/tex" id="MathJax-Element-91">O(n^2)</script>。

动态加点?
每次增加一个点,一般求法会重新计算一次达到 O(n2) O ( n 2 ) <script type="math/tex" id="MathJax-Element-92">O(n^2)</script>的复杂度,而重心拉格朗日插值法只需 O(n) O ( n ) <script type="math/tex" id="MathJax-Element-93">O(n)</script>计算 wi w i <script type="math/tex" id="MathJax-Element-94">w_i</script>即可。

横坐标连续?
不需多次计算 wi w i <script type="math/tex" id="MathJax-Element-95">w_i</script>,单次求值 O(n) O ( n ) <script type="math/tex" id="MathJax-Element-96">O(n)</script>,预处理逆元的话求原函数也只会 O(n) O ( n ) <script type="math/tex" id="MathJax-Element-97">O(n)</script>。

应用

bzoj4559:成绩比较


xi=1UxiRxUxinRx1(Ux109,Rx102,n102) ∏ x ∑ i = 1 U x i R x ( U x − i ) n − R x − 1 ( U x ≤ 10 9 , R x ≤ 10 2 , n ≤ 10 2 )
<script type="math/tex; mode=display" id="MathJax-Element-43">\prod_x\sum_{i=1}^{U_x}i^{R_x}(U_x-i)^{n-R_x-1}(U_x\le 10^9,R_x\le 10^2,n\le 10^2)</script>

很显然是一个 n n <script type="math/tex" id="MathJax-Element-44">n</script>次多项式,对于每一个 x x <script type="math/tex" id="MathJax-Element-45">x</script>计算 n+1 n + 1 <script type="math/tex" id="MathJax-Element-46">n+1</script>个点值确定函数,然后带入 Ux U x <script type="math/tex" id="MathJax-Element-47">U_x</script>计算即可,因为 xi x i <script type="math/tex" id="MathJax-Element-48">x_i</script>连续,所以时间复杂度 O(n2) O ( n 2 ) <script type="math/tex" id="MathJax-Element-49">O(n^2)</script>。

bzoj2655: calc

听说结果是次数为 2n 2 n <script type="math/tex" id="MathJax-Element-50">2n</script>的多项式,那么先DP出小的点,直接上拉格朗日就行了。

bzoj3453:XLkxc

f(n)=i=0nj=1a+ibx=1jxk(k123,n,a,b109) f ( n ) = ∑ i = 0 n ∑ j = 1 a + i b ∑ x = 1 j x k ( k ≤ 123 , n , a , b ≤ 10 9 )
<script type="math/tex; mode=display" id="MathJax-Element-73">f(n)=\sum_{i=0}^{n}\sum_{j=1}^{a+ib}\sum_{x=1}^{j}x^k(k\le 123,n,a,b\le 10^9)</script>

虽然看上去很复杂,其实拆开并不复杂。
首先:

x=1jxk ∑ x = 1 j x k
<script type="math/tex; mode=display" id="MathJax-Element-74">\sum\limits_{x=1}^{j}x^k</script>是一个 k+1 k + 1 <script type="math/tex" id="MathJax-Element-75">k+1</script>次多项式(二项式定理展开后差分)。

其次:

g(m)=j=1mx=1jxk g ( m ) = ∑ j = 1 m ∑ x = 1 j x k
<script type="math/tex; mode=display" id="MathJax-Element-76">g(m)=\sum_{j=1}^{m}\sum_{x=1}^jx^k</script>
是一个 k+2 k + 2 <script type="math/tex" id="MathJax-Element-77">k+2</script>次多项式(同理)。

可以对 g(m) g ( m ) <script type="math/tex" id="MathJax-Element-78">g(m)</script>进行插值求大的 m m <script type="math/tex" id="MathJax-Element-79">m</script>数。

最后:

f(n)=i=0ng(a+ib) f ( n ) = ∑ i = 0 n g ( a + i b )
<script type="math/tex; mode=display" id="MathJax-Element-80">f(n)=\sum_{i=0}^ng(a+ib)</script>是一个 k+3 k + 3 <script type="math/tex" id="MathJax-Element-81">k+3</script>次多项式(同理)。
再对 f(n) f ( n ) <script type="math/tex" id="MathJax-Element-82">f(n)</script>进行插值就好了。

而且得出了两个结论:
1.每个求和号一般会增加多项式次数1.
2.插值可以嵌套。
特别是第二点,基本上可以在 O(k2) O ( k 2 ) <script type="math/tex" id="MathJax-Element-83">O(k^2)</script>解决所有跟多项式有关的题。

Code for bzoj3453:

#include<bits/stdc++.h>
using namespace std;
struct IO
{
    streambuf *ib,*ob;
    int buf[50];
    inline void init()
    {
        ios::sync_with_stdio(false);
        cin.tie(NULL);cout.tie(NULL);
        ib=cin.rdbuf();ob=cout.rdbuf();
    }
    inline int read()
    {
        char ch=ib->sbumpc();int i=0,f=1;
        while(!isdigit(ch)){if(ch=='-')f=-1;ch=ib->sbumpc();}
        while(isdigit(ch)){i=(i<<1)+(i<<3)+ch-'0';ch=ib->sbumpc();}
        return i*f;
    }
    inline void W(int x)
    {
        if(!x){ob->sputc('0');return;}
        if(x<0){ob->sputc('-');x=-x;}
        while(x){buf[++buf[0]]=x%10,x/=10;}
        while(buf[0])ob->sputc(buf[buf[0]--]+'0');
    }
}io;
const int mod=1234567891,lim=130;
int k,a,n,d;
int fg[150],fy[150],wf[150],fac[150],inv[150];
inline int power(int a,int b)
{
    int res=1;
    for(;b;b>>=1)
    {
        if(b&1)res=1ll*res*a%mod;
        a=1ll*a*a%mod;
    }
    return res;
}
int DIV;
inline int calc(int *w,int u)
{
    if(u<=lim)return w[u];
    static int l,res,Div;
    l=1,res=0,Div=DIV;
    for(int i=1;i<=lim;i++)l=1ll*l*((1ll*u-i+mod)%mod)%mod;
    for(int i=1;i<=lim;i++)
    {
        res=(1ll*res+1ll*Div*power(((1ll*u-i+mod)%mod),mod-2)%mod*w[i]%mod)%mod;
        Div=1ll*Div*(i-lim+mod)%mod*inv[i]%mod;
    }
    return 1ll*res*l%mod;
}
int main()
{
    io.init();int T=io.read();
    fac[0]=fac[1]=inv[0]=inv[1]=1;DIV=1;
    for(int i=2;i<=lim;i++)inv[i]=(-1ll*(mod/i)*inv[mod%i]%mod+mod)%mod;
    for(int i=1;i<=lim;i++)fac[i]=1ll*fac[i-1]*i%mod;
    for(int i=2;i<=lim;i++)DIV=1ll*DIV*power(1-i+mod,mod-2)%mod;
    while(T--)
    {
        k=io.read(),a=io.read(),n=io.read(),d=io.read();
        for(int i=1;i<=lim;i++){fg[i]=power(i,k);}
        for(int i=1;i<=lim;i++){fg[i]=(1ll*fg[i]+fg[i-1])%mod;}
        for(int i=1;i<=lim;i++){fg[i]=(1ll*fg[i]+fg[i-1])%mod;}
        fy[0]=calc(fg,a);
        for(int i=1;i<=lim;i++)
        {
            fy[i]=(1ll*fy[i-1]+calc(fg,(1ll*a+1ll*i*d)%mod))%mod;
        }
        printf("%d\n",calc(fy,n));
    }
}
©️2020 CSDN 皮肤主题: 编程工作室 设计师:CSDN官方博客 返回首页