1. 程式人生 > 實用技巧 >lagrange 插值法模板

lagrange 插值法模板

對於n次多項式

\[\Gamma(k) = \sum_{i=0}^{n}y_i\prod_{j=0,{j}\neq {i} }^{n}\frac{k-x_j}{x_i-x_j} \]

const int mod=1e9+7;
template<class T>
T qpow(T x,int k){
    T e=1;
    while(k){
        if(k&1)e=1ll*e*x%mod;
        x=1ll*x*x%mod;
        k>>=1;
    }
    return e;
}
template<class T>
T lagrange(int n,int x[],int y[],T k){
    T ans=0;
    for(int i=0;i<=n;++i){
        T fz=1,fm=1;//分子分母
        for(int j=0;j<=n;++j){
            if(i==j)continue;
            fm=1ll*fm*(x[i]-x[j])%mod;
            fz=1ll*fz*(k-x[j])%mod;
        }
        ans=(1ll*ans+1ll*y[i]*fz%mod*qpow(fm,mod-2)%mod)%mod;
    }
    return (ans+mod)%mod;
}

在x取值連續時的做法

\(x_i\)​ 的取值都是連續的,可以把上面的演算法優化到O(n)複雜度
首先把\(x_i\)換成i,新的式子為

\[f(k) = \sum_{i=0}^{n} y_i \prod_{i \not = j} \frac{k - j}{i - j}f(k) \]

對於 \(\prod_{i \not = j} ^{n}\frac{k - j}{i - j}​\) 可以優化
對於分子來說,我們維護出關於k的字首積和字尾積,也就是

\[pre_i = \prod_{j = 0}^{i-1} k - j \]

\[suf_i = \prod_{j = i+1}^n k - j \]

對於分母來說,觀察發現這其實就是階乘的形式,我們用fac[i]來表示 \(i!\)
那麼式子就變成

\[f(k) = \sum_{i=0}^n y_i \frac{pre_i * suf_i }{fac[i] * fac[N - i]}$$​ 分母可能會出現符號問題,也就是說,當N - i為奇數時,分母應該取負號 ```cpp const int mod=1e9+7; template<class T> T qpow(T x,int k){ T e=1; while(k){ if(k&1)e=1ll*e*x%mod; x=1ll*x*x%mod; k>>=1; } return e; } int lagrange(int n,int x[],int y[],int k){ int pre[n+1],suf[n+1],fac[n+1],ans=0; fac[0]=pre[0]=suf[n]=1; for(int i=1; i<=n; ++i)pre[i]=1ll*(k-x[i-1])*pre[i-1]%mod; for(int i=n-1;i>=0;++i)suf[i]=1ll*(k-x[i+1])*suf[i+1]%mod; for(int i=1;i<=n;++i)fac[i]=1ll*i*fac[i-1]%mod; for(int i=0;i<=n;++i){ int fz=1ll*pre[i]*suf[i]%mod;//分子分母 int fm=1ll*fac[i]*fac[n-i]%mod; if((n-i)&1)fm=(mod-fm)%mod; ans=(1ll*ans+1ll*y[i]*fz%mod*qpow(fm,mod-2)%mod)%mod; } return (ans+mod)%mod; } ``` ### 2019南昌邀請賽B. Polynomial [題目連結](https://nanti.jisuanke.com/t/40254) Define $f_{(x)} = a_0+a_1 \times x^1 + a_2 \times x^2 +…+ a_n \times x^n$ the number maybe tremendous huge, so take the module of 9999991. For a polynomial, XH doesn't know any $a_i$ ​ , but he knows the value of $f_{(0)} , f_{(1)}, f_{(2)}…f_{(n)}$ . He wants to calculate $\displaystyle\sum_{i=L}^R f_{(i)}$ mod 9999991. #### Input The first line contains a number $T ( 1 \le T \le 5 )$ ——the number of test cases. For each test case. The first line contains two integers $n (1 \le n \le 1000)$ and $m(1 \le m \le 2000)$. The second line contains n+1 integers,representing the value of $f_{(0)} , f_{(1)}, f_{(2)}…f_{(n)}$ ​ The next mm lines, each line contains two integers L, R. $( 1 \le L \le R \le 9999990)$ #### Output For each test case, output m lines, each line contains one integer, the value of $isplaystyle\sum_{i=L}^R f_{(i)}$ mod 9999991. #### 樣例輸入 ``` 1 3 2 1 10 49 142 6 7 95000 100000 ``` #### 樣例輸出 ``` 2519 1895570 ``` 題意就是知道(0, $f_0$) , (1, $f_1$) ... (n, $f_n$) n+1個點,求$S(k)=\sum_{i=0}^k f(i)$ 由n+1個點可求得 $f(n+1)$ 再由$S(k)-S(k-1)=f(k)$可知,$S(k)=\sum_{i=0}^k f(i)$是n+1次多項式 且我們可以求出這n+2個點,分別是 (0, $f_0$) , (1, $f_0+f_1$), ..., (n, $f_0+f_0+... +f_n$), (n, $f_0+f_0+... +f_n+f_{n+1})$ 這樣,求出$S(k)$ 之後 $S(r)-S(l-1)$即為答案。 程式碼: ```cpp #include<bits/stdc++.h> using namespace std; char buf[100000],*L=buf,*R=buf; #define gc() L==R&&(R=(L=buf)+fread(buf,1,100000,stdin),L==R)?EOF:*L++ template<typename T> inline void read(T&x) { int flag=x=0; char ch=gc(); while (ch<'0'||ch>'9')flag|=ch=='-',ch=gc(); while (ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+ch-48,ch=gc(); if(flag)x=-x; } #define Init(arr,val) memset(arr,val,sizeof(arr)) const int inf=0x3f3f3f3f,mod=9999991,MAXN=1020; typedef long long ll;//鬼畜資料,只能全部用ll,不然一個點都過不了,浪費了了我一下午。。。。 template<class T> T qpow(T x,int k) { T e=1; while(k) { if(k&1)e=1ll*e*x%mod; x=1ll*x*x%mod; k>>=1; } return e%mod; } ll fac[MAXN],ffa[MAXN]; //fac是階乘 //ffa是分母,對於m次詢問,都是相同的分母,儲存起來, //因為S(x) 已經確定,分母裡的快速冪只需求一次即可,不然坑爹資料一個點也過不了 ll lagrange(ll n,ll x[],ll y[],ll k) { ll pre[MAXN],suf[MAXN],ans=0; fac[0]=pre[0]=suf[n]=1; for(int i=1; i<=n; ++i)pre[i]=1ll*(k-x[i-1])*pre[i-1]%mod; for(int i=n-1; i>=0; --i)suf[i]=1ll*(k-x[i+1])*suf[i+1]%mod; for(int i=0; i<=n; ++i) { ll fz=1ll*pre[i]*suf[i]%mod; ans=(1ll*ans+1ll*y[i]*fz%mod*ffa[i]%mod)%mod; } return (ans+mod)%mod; } int main() { fac[0]=1; for(int i=1; i<MAXN; ++i)fac[i]=1ll*i*fac[i-1]%mod; ll n,m,x[MAXN],y[MAXN]; int t; read(t); while(t--) { read(n),read(m); for(int i=0; i<=n; ++i) { x[i]=i; read(y[i]); } for(int i=0; i<=n; ++i) {//f(x) ll fm=1ll*fac[i]*fac[n-i]%mod; if((n-i)&1)fm=(mod-fm)%mod; ffa[i]=1ll*qpow(fm,mod-2)%mod; } y[n+1]=lagrange(n,x,y,n+1);//f(n+1) x[n+1]=n+1; for(int i=1; i<=n+1; ++i)y[i]+=y[i-1];//字首和就是S(x)個點的y值 for(int i=0; i<=n+1; ++i) {//m次查詢,分母記錄下來,求一次即可 ll fm=1ll*fac[i]*fac[n+1-i]%mod; if((n+1-i)&1)fm=(mod-fm)%mod; ffa[i]=1ll*qpow(fm,mod-2)%mod; } while(m--) { ll l,r; read(l),read(r); printf("%lld\n",(lagrange(n+1,x,y,r)-lagrange(n+1,x,y,l-1)+mod)%mod); } } return 0; } ```\]