拉格朗日插值模板(P4781)
阿新 • • 發佈:2021-06-28
題目連結:拉格朗日插值
拉格朗日插值:給定 \(k+1\) 個點對 \((x_i,y_i)\) (\(x_i\)各不相同)能夠唯一確定一個最高次為 \(k\) 次的多項式,那麼如何進行構造,來求該多項式呢?我們先以經過 \((x_1,1),(x_2,0),(x_3,0)\) 這三個點的4次多項式為例:那麼我們可以進行構造設 \(f(X)=\frac{(X-x_2)*(X-x_3)}{(x_1-x_2)*(x_1-x_3)}\),我們帶入這三個點發現顯然該多項式是我們所求的多項式。那麼如果我們將該三個點擴充套件到 \((x_1,y_1),(x_2,y_2),(x_3,y_3)\) 呢?
顯然:
- \(f_1(X)=\frac{(X-x_2)*(X-x_3)}{(x_1-x_2)*(x_1-x_3)}*y_1\)
- \(f_2(X)=\frac{(X-x_1)*(X-x_3)}{(x_2-x_1)*(x_2-x_3)}*y_2\)經過\((x_2,y_2)\),且\(f_2(x_1) = 0,f_2(x_3) = 0\)
- \(f_3(X)=\frac{(X-x_1)*(X-x_2)}{(x_3-x_1)*(x_3-x_2)}*y_3\)經過\((x_3,y_3)\),且\(f_3(x_1) = 0,f_3(x_2) = 0\)
所以我們最後構造的多項式為\(F(X) = f_1(X) + f_2(X) + f_3(X)\)
這樣就完成了拉格朗日插值;模板:
\(Code:\)
/* -*- encoding: utf-8 -*- ''' @File : 拉格朗日模板.cpp @Time : 2021/06/28 11:21:28 @Author : puddle_jumper @Version : 1.0 @Contact : [email protected] ''' # here put the import lib*/ #include<set> #include<iostream> #include<cstring> #include<cmath> #include<cstdio> #include<cstdlib> #include<map> #include<algorithm> #include<vector> #include<queue> #define ch() getchar() #define pc(x) putchar(x) #include<stack> #include<unordered_map> #define rep(i,a,b) for(auto i=a;i<=b;++i) #define bep(i,a,b) for(auto i=a;i>=b;--i) #define lowbit(x) x&(-x) #define ll long long #define ull unsigned long long #define pb push_back #define mp make_pair #define PI acos(-1) using namespace std; template<typename T>void read(T&x){ static char c; static int f; for(c=ch(),f=1; c<'0'||c>'9'; c=ch())if(c=='-')f=-f; for(x=0; c>='0'&&c<='9'; c=ch())x=x*10+(c&15); x*=f; } template<typename T>void write(T x){ static char q[65]; int cnt=0; if(x<0)pc('-'),x=-x; q[++cnt]=x%10,x/=10; while(x) q[++cnt]=x%10,x/=10; while(cnt)pc(q[cnt--]+'0'); } const int N = 2e3+10; const ll mod = 998244353; int n;ll k; ll x[N],y[N]; ll qpow(ll a,ll b){ ll res = 1ll; while(b){ if(b&1){ res *= a;res%=mod; } b>>=1; a*=a;a%=mod; } return res%mod; } ll inv(ll x){ return qpow(x,mod-2); } void solve(){ read(n);read(k); rep(i,1,n)read(x[i]),read(y[i]); ll ans = 0ll; rep(i,1,n){ ll res = y[i]; rep(j,1,n){ if(i == j)continue; res *= (k-x[j]); res %= mod; res *= inv(x[i] - x[j]); res %= mod; } ans += (res%mod + mod)%mod; ans %= mod; } write(ans);pc('\n'); } signed main(){solve();return 0; }