一点点拉格朗日插值

抄起我的混凝土数学，复习了会上升下降积分，看了看斯特林数和伯努利数。结果误把伯努利数的递推式看成了O(n)O(n)的，改了半天一到3次方以上就翻车，这才反应过来。

8会啊…找了很长时间才知道是拉格朗日插值。

有这么一条定理。

平面上的n+1个点可以确定一个n次多项式。

而∑ip\sum i^p是一个p+1p+1次多项式，况且在p≤1000000p \leq 1000000的前提下很容易找到前p+1p+1个点。如果我们找到一种方式去构造这样一个多项式问题就解决了。

拉格朗提出了这么一种方法。

构造

设函数

ℓj(x)=∏k≠jx−xkxj−xk\ell_j(x)=\prod_{k \neq j} \frac{x-x_k}{x_j-x_k}

注意其特性

  • xjx_j被带入时，其为11
  • 其他任何xx被带入时，会出现00项。

因此这个东西就可以为任何一个点对创建在多项式中的单元yjℓj(xj)y_j \ell_j(x_j)。剩下的问题是如何快速求解。

当我们只求单点（率先带入点），且整数x连续取值时，观察到分母和分子都变为阶乘形式，且分式的构成元素在前后两项之间具有递推的关系。

ℓj(x)=(x−1)(x−2)⋯(x−(j−1))(x−(j+1))⋯(x−n)(j−1)(j−2)⋯(j−(j−1))(j−(j+1))⋯(j−n))\ell_j(x)=\frac{(x-1)(x-2)\cdots (x-(j-1))(x-(j+1))\cdots (x-n)}{(j-1)(j-2) \cdots (j-(j-1))(j-(j+1)) \cdots (j-n))}

可以看到，分子和分母项都能够拆分为左右两个部分。分别维护前缀后缀与阶乘逆元即可。

阶乘并不需要全部的，仅仅是左侧和右侧的一小部分。

板子

嫖的，自己的写爆了而且巨丑。

1
2
3
4
5
6
7
8
9
10
11
12
int lagrange(int n,int *x,int *y,int k) {
    int ans=0;
    for(int i=1;i<=n;++i) {
        int f=1,g=1;
        for(int j=1;j<=n;++j) if(i!=j) {
            f=1LL*f*(k-x[j]+mod)%mod;
            g=1LL*g*(x[i]-x[j]+mod)%mod;
        }
        upd(ans,1LL*y[i]*f%mod*inv(g)%mod);
    }
    return ans;
}
1
2
3
4
5
6
7
8
9
10
11
12
int lagrange(int n,int *y,int k) {
    int ans=0;
    pre[0]=suf[n+1]=1;
    for(int i=1;i<=n;++i) pre[i]=1LL*pre[i-1]*(k-i)%mod;
    for(int i=n;i>=1;--i) suf[i]=1LL*suf[i+1]*(k-i)%mod;
    for(int i=1;i<=n;++i) {
        int a=1LL*pre[i-1]*suf[i+1]%mod*inv[i-1]%mod*inv[n-i]%mod;
        if((n-i)&1) a=mod-a;
        upd(ans,1LL*a*y[i]%mod);
    }
    return ans;
}

错误的代码

这份看错了伯努利数的递推复杂度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#include <cstdio>
#include <algorithm>
#include <string>
#include <vector>
#include <cstring>
#include <iostream>
#include <queue>
#include <set>
#include <cmath>
#include <map>
using namespace std;
using ll=long long;
const int MAXN=1000010;
const int mod=1e9+7;


ll inv[MAXN];
double B[MAXN];
double C[MAXN];

void init(int n){
    //niyuan
    inv[1]=1;
    for(int i=2;i<=n+1;i++) {
        inv[i] = inv[mod % i] * (mod - mod / i) % mod;
    }

    //C
    C[0]=1;
    for(int i=1;i<=n+1;i++){
//        C[i]=C[i-1]*(n+1-i+1)%mod*inv[i]%mod;
        C[i]=C[i-1]*(n+1-i+1)/i;
    }
    
    
    //WRONG
    //伯努利数的递推是n^2级别的，理解式子上出现了失误。
    B[0]=1;
    double s=C[0];
    for(int i=1;i<=n;i++){
//        B[i]=-inv[i+1]*s%mod;
//        s=(s+B[i]*C[i]%mod)%mod;
        B[i]=-s/(i+1);
        s=(s+B[i]*C[i]);
    }
}

ll n2[MAXN];
int main(){
    ll n,m;cin>>n>>m;

    init(m);
    for(int i=0;i<=m;i++)cout<<C[i]<<" ";
    cout<<endl;
    for(int i=0;i<=m;i++)cout<<B[i]<<" ";
    cout<<endl;

    n2[0]=1;
    for(int i=1;i<=m+1;i++){
        n2[i]=n2[i-1]*n%mod;
    }

    for(int i=0;i<=m+1;i++)cout<<n2[i]<<" ";
    cout<<endl;

    double ans=0;
    for(int i=0;i<=m;i++){
//        ans=(ans+C[i]*B[i]%mod*n2[m+1-i]%mod)%mod;
        ans=(ans+C[i]*B[i]*n2[m+1-i]);
    }
//    (ans*=inv[m+1])%=mod;
    ans=ans/(m+1);
    cout<<ans<<endl;


}

这份自己yy的拉格朗日插值在阶乘的求法上有问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include <cstdio>
#include <algorithm>
#include <string>
#include <vector>
#include <cstring>
#include <iostream>
#include <queue>
#include <set>
#include <cmath>
#include <map>
using namespace std;
using ll=long long;
const int MAXN=1000010;
const int mod=1e9+7;

ll fac[MAXN];
ll facinv[MAXN];
ll inv[MAXN];

ll qpow(ll a,ll b,ll p=mod){
    ll res=1;
    for(;b;b>>=1,a=a*a%p){
        if(b&1){
            res=res*a%p;
        }
    }
    return res;
}

inline ll get_inv(ll a,ll p=mod){
    return qpow(a,p-2,p);
}

inline ll Fac(int n,bool pos=true){
    int flg=1;
    if(!pos && n%2==1)flg=-1;
    return fac[n]*flg;
}

inline ll Invfac(int n,bool pos=true){
    int flg=1;
    if(!pos && n%2==1)flg=-1;
    return facinv[n]*flg;
}


void init(int n){
    fac[0]=1;
    for(int i=1;i<=n;i++)fac[i]=fac[i-1]*i%mod;
    facinv[n]=get_inv(fac[n]);
    for(int i=n-1;i>=1;i--){
        facinv[i]=facinv[i+1]*(i+1)%mod;
    }
    facinv[0]=1;

    inv[1]=1;
    for(int i=2;i<=n;i++)
        inv[i] = inv[mod % i] * (mod-mod/i) % mod;
}

ll ans=0;
ll yy[MAXN];

int main(){
    int n,m;cin>>n>>m;

    int M=m+2;
    init(max(n,M));
    for(int i=1;i<=M;i++){
        yy[i]=qpow(i,m);
    }
    for(int i=1;i<=M;i++){
        yy[i]=(yy[i]+yy[i-1])%mod;
    }
    for(int i=1;i<=M;i++)cout<<yy[i]<<" ";
    cout<<endl;

    if(n<=M){
        cout<<yy[n]<<endl;
        return 0;
    }


    for(int i=1;i<=M;i++){
        ll cell=yy[i]*Fac(n-1)%mod*Invfac(n-M-1)%mod*inv[n-i]%mod;
        (cell*=Invfac(i-1))%=mod;
        (cell*=Invfac(abs(i-M),false))%=mod;
        (ans+=cell)%=mod;
//        ll cell=yy[i]*Fac(n-1)/Fac(n-M-1)/(n-i);
//        cell/=Fac(i-1);
//        cell/=Fac(abs(i-M),false);
//        ans+=cell;
    }

    cout<<(ans+mod)%mod<<endl;
}
学习
一点点拉格朗日插值
https://blog.chenc.me/2020/03/19/lagrange-polynomial/
作者
CC
发布于
2020年3月19日
许可协议