拉格朗日插值

前言 : 本人数学水平较菜, 若有出错, 稍加原谅。

简介 :

一种由 \(n + 1\) 个点 \((x_i, y_i)\) 去求出最多为 \(n\) 次的多项式 : $y = f(x) $ 使得其经过每一个点, 或者求出另一个在此函数的点 \((x, y)\)

特定的, 我们钦定下文的 \(i, j \in[0, n]\)

拉格朗日基本多项式

\[l_i(x) = \prod_{j\not=i}^n \frac{x-x_i}{x_i-x_j} \]

因为其构造巧妙, 我们可以发现 \(l_i(x_i) = 1\), \(l_i(x_j) = 0\) 这一奇妙性质。

那么继续构造 \(f(x)\) 可以得到 :

\[f(x) = \sum_{i=0}^ny_il_i(x) \]

所以我们可以在 \(O(n)\) 的时间内计算单个 \(l_i(x)\) , 从而实现 \(O(n^2)\) 插值出新的点 \((x, y)\)

重心拉格朗日插值

这是一种增量式的拉格朗日插值计算方法。

考虑展开 \(f(x)\) :

\[f(x) = \sum_{i=0}^ny_i\prod_{j=0, j\not=i}^n \frac{x-x_i}{x_i-x_j} \]

我们令 \(g =\prod_{i=0}^n x-x_i\) ,原式为 :

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

我们只需要计算 \(t_i = \prod_{i\not=j}^n\frac{y_i}{(x-x_i)(x_i-x_j)}\) 就可以得到 :

\[f(x) =g\sum_{i = 0} ^ n \frac{t_i}{x - x_i} \]

我们发现 \(t_i\) 是可以 \(O(n)\) 计算的, 十分优秀。

\(\rm Bonus :\) 取值连续

很多时候, \(x_i\in[1, n]\) 这段区间, 我们发现拉格朗日插值可以继续优化。

\[f(x) = \sum_{i=0}^ny_i\prod_{j\not=i}^n \frac{x-x_i}{x_i-x_j} \]

这个式子, 我们记录 :

\[pre_i = \prod_{j=0}^nx-x_j \]

\[suf_i = \prod_{j=i}^nx-x_j \]

对此, 我们发现非常巧妙的, 分母是阶乘的形式 :

\[f(x)=\sum_{i=0}^{n}y_i\frac{pre_{i-1}*suf_{i + 1}}{fac_{i -1}*fac_{n - i}} \]

注意, 若 \(n -i\) 位奇数, 那么分母是负的。

自此拉格朗日插值可以优化到 \(O(n\log n + nt)\) 其中 求逆元 + 单点求值。

板子 :

namespace lagrange {
	
	// nNode -> n - 1 -> f(x)
	const int MOD = 1e9 + 7, N = 5e3 + 5;
	
	int n, k, x[N], y[N], pre[N], suf[N], fac[N];
	
	inline void mul (int & x, int y) { x = x * y % MOD; }
	inline void add (int & x, int y) { x = (x + y) % MOD; }
	inline int qpow (int x, int y) {
		int s = 1;
		for ( ; y; mul(x, x), y >>= 1)
			if (y & 1) mul(s, x);
		return s;
	}
	
	inline int getG (int i, int pos) {
		int ans = y[i], s1 = 1, s2 = 1;
		for (int j = 1; j <= n; j ++) 
			if (j != i) mul(s1, pos - x[j]), mul(s2, x[i] - x[j]);
		ans = ans * (s1 * qpow(s2, MOD - 2) % MOD) % MOD;
		return ans;
	}
	inline void extendGet () {
		int ans = 0;
		cin >> n >> k;
		for (int i = 1; i <= n; i ++) cin >> x[i] >> y[i];
		for (int i = 1; i <= n; i ++) add(ans, getG(i, k));
		ans = (ans + MOD) % MOD, cout << ans;
	}
	inline int simpleGet (int n, int k) {
		// k + 1
		int ans = 0, y = 0;
		pre[0] = suf[k + 3] = fac[0] = 1;
		for (int i = 1; i <= k + 2; i ++) pre[i] = pre[i - 1] * (n - i) % MOD;
		for (int i = k + 2; i >= 1; i --) suf[i] = suf[i + 1] * (n - i) % MOD;
		for (int i = 1; i <= k + 2; i ++) fac[i] = fac[i - 1] * i % MOD;
		for (int i = 1; i <= k + 2; i ++) {
			add(y, qpow(i, k)); // get pos-val
			int a = pre[i - 1] * suf[i + 1] % MOD, b = fac[i - 1] * ((k - i) & 1 ? -1 : 1) * fac[k + 2 - i] % MOD;
			add(ans, y * a % MOD * qpow(b, MOD - 2) % MOD);
		}
		ans = (ans + MOD) % MOD;
		return ans;
	}
}

习题 :

P4781 模板

CF622F

\(\sum_i^n i^k\) 是个 \(k + 1\) 次多项式, 那我们拿 \([1, k+ 2]\) 去插值就可以了。

P4593

\(\sum\sum_i i^{m+1}\) 大力拉插。

P3270

Orz 牛逼组合题。

我们记 \(f_{i, j}\) 为前 \(i\) 门碾压 \(j\) 个人, 考虑到 \(j\) 单调不升, 有如下转移 :

\[f_{i, j} = \sum_{t=j}^kf_{i-1,t}\tbinom{t}{t-j}\tbinom{n-k-1}{r_i - 1 - (k - j)}g(i) \]

其中 \(g(i)\) 为这一门课的方案数, 枚举划分即可。

\[g(i) = \sum_{j = 1}^{u_i}j^{n-r_i}(u_i-j)^{r_i-1} \]

其中 \(g(i)\) 可以直接插值。

\[\rm O(n^2m \log_2 u_{max}) \]

解决这一问题。

  • 差分理论 :

\(f(x)\)\(n\) 次的, 那么 \(f(x) - f(x - 1)\)\(n - 1\) 次的。

P4463 :

\[f_{i, j} - f_{i, j - 1} = j * f_{i - 1, j - 1} \]

猜测为关于 \(j\) 的多项式 :

\[g(i) - 1= g(i - 1) + 1 \]

\[g(i) = g(i - 1) + 2 \]

那么就可以插值了。

CF995F ;

\[f_{u, d} = f_{u, d - 1} + \prod f_{v, d} \]

\[f_{u, d} - f_{u, d - 1} = \prod f_{v, d} \]

考虑这是关于 \(d\) 的函数 :

\[g_u - 1= \sum g_v \]

我不会归纳之类的, 但是考虑到叶子结点 \(g_v = 1\)

所以对于其父亲 \(u\) 有 : \(g_u = size_u\)

那么这样一层一层地推上去,都有 : \(g_{rt} = size_{rt}\)

所以, \(f_{1, d}\) 是一个关于 \(d\)\(size_1\) 次多项式。

大力插值即可。

posted @ 2022-07-12 00:00  Cust10  阅读(176)  评论(0)    收藏  举报