• 【笔记】拉格朗日插值还原系数表达式


    大家都知道拉格朗日插值的公式(已知 \(n\) 个点 \((x_i,y_i)\),求唯一确定的经过这 \(n\) 个点的\(n-1\) 次多项式):

    \[f(x)=\sum_{i=1}^ny_i\left(\prod_{j\ne i}\frac{x-x_j}{x_i-x_j}\right) \]

    但是洛谷模板只让求了一个点的点值,就不用把系数表达式算出来。本蒻稽不会 \(O(n\log n)\) 的多项式快速插值又想在省选骗分,只好学一下 \(O(n^2)\) 的插值方法。

    上面公式分成两部分,一部分是分子 \(y_i \prod_{j\ne i} \frac 1{x_i-x_j}\)\(O(n)\) 可以计算,但分母部分 \(\prod_{j\ne i}(x-x_j)\) 需要 \(O(n^2)\)

    优化方案是先把 \(F(x)=\prod_{i=1}^n (x-x_i)\) 算出来,然后用 \(F(x)/(x-x_j)\) 来得到 \(f_j(x)\).

    \(F(x)_i\)\(F(x)\)\(i\) 项的系数,\(f_j(x)_i\) 同理,有:

    \[F(x)_0=-x_jf_j(x)_0 \]

    \[F(x)_i=f_j(x)_{i-1}-x_jf_j(x)_{i},\quad (i>1) \]

    于是得到:

    \[f_j(x)_0=-F(x)_0/x_j \]

    \[f_j(x)_i=(f_j(x)_{i-1}-F(x)_i)/x_j,\quad (i>1) \]

    但是注意这样手动模拟的前提是 \(x_j\ne 0\)\(x_j=0\) 的情况要特判掉(\(f_j(x)=F(x)/x\)).

    算出 \(f_j(x)\) 后,把之前的分子部分乘到系数上。最后 \(f(x)=\sum_{j=1}^n f_j(x)\) 就做完了。

    然后把代码写出来就好了(洛谷模板):

    #include <bits/stdc++.h>
    
    using namespace std;
    typedef long long ll;
    const int N = 2005, D = 998244353;
    int fpm(int a, int b = D - 2, int c = D) {
    	int res = 1 % c;
    	while(b) {
    		if(b & 1) res = (ll)res * a % c;
    		b >>= 1, a = (ll)a * a % c;
    	}
    	return res;
    }
    inline int _(int a) { return a + (a >> 31 & D); }
    inline void Add(int &a, ll b) { a = (a + b) % D; }
    inline void Sub(int &a, ll b) { a = _((a - b) % D); }
    inline void Div(int &a, int b) { a = (ll)a * fpm(b) % D; }
    inline void Mul(int &a, int b) { a = (ll)a * b % D; }
    int n, k, x[N], y[N], f[N];
    void trans(int x[], int y[], int f[]) {
    	static int c[N], F[N], t[N], ix[N];
    	for(int i = 1; i <= n; i++) c[i] = 1, f[i] = 0;
    	for(int i = 1; i <= n; i++) {
    		ix[i] = fpm(x[i]);
    		for(int j = 1; j <= n; j++)
    			if(j != i) Mul(c[i], _(x[i] - x[j]));
    		c[i] = (ll)fpm(_(c[i])) * y[i] % D;
    	}
    	for(int i = 0; i <= n; i++) F[i] = 0;
    	F[0] = 1;
    	for(int i = 1; i <= n; i++)
    		for(int j = i; j >= 0; j--) 
    			F[j] = _((F[j - 1] - (ll)x[i] * F[j]) % D);
    	for(int i = 1; i <= n; i++) {
    		if(!x[i]) for(int j = 1; j <= n; j++) t[j - 1] = F[j];
    		else {
    			t[0] = _((ll)-F[0] * ix[i] % D);
    			for(int j = 1; j < n; j++)
    				t[j] = _((ll)(t[j - 1] - F[j]) * ix[i] % D);
    		}
    		for(int j = 0; j < n; j++) Add(f[j], (ll)t[j] * c[i] % D);
    	}
    }
    
    int main() {
    	ios::sync_with_stdio(false); cin.tie(nullptr);
    	cin >> n >> k;
    	for(int i = 1; i <= n; i++) cin >> x[i] >> y[i];
    	trans(x, y, f);
    	int tmp = 1, ans = f[0];
    	for(int i = 1; i < n; i++)
    		Mul(tmp, k), Add(ans, (ll)tmp * f[i]);
    	cout << ans << '\n';
    
    	return 0;
    }
    
    
  • 相关阅读:
    洛谷P3959 宝藏(状压dp)
    洛谷P3645 [APIO2015]雅加达的摩天楼(最短路+分块)
    洛谷P3646 [APIO2015]巴厘岛的雕塑(数位dp)
    洛谷P4770 [NOI2018]你的名字(后缀自动机+线段树)
    洛谷P4768 [NOI2018]归程(克鲁斯卡尔重构树+最短路)
    hive3.1.1 hive-site.xml
    mysql 远程连接数据库的二种方法
    linux彻底干干净净完全卸载 mysql
    jdk环境变量配置
    Eclipse常用快捷键
  • 原文地址:https://www.cnblogs.com/huaruoji/p/15926926.html
Copyright © 2020-2023  润新知