闲话 23.1.7

闲话

数学这里证明是个避不过去的东西。你总得面对他。

—— 巨佬

今日 ABC（
先写一点 慢慢补上

upd: 挺耻辱赛的感觉

拉格朗日插值

熟知平面上 \(n + 1\) 个点可以唯一确定一个 \(n\) 次多项式。可是这系数要怎么确定呢？我们已经有了高斯消元法，可是其复杂度是 \(O(n^3)\) 的，不优秀。拉格朗日提供了一种复杂度为 \(O(n^2)\) 的做法，其思想就是直接配凑。

考虑平面上的 \(n\) 个点，第 \(i\) 个点记作 \((x_i, y_i)\)，并记 \(F(x)\) 为这 \(n\) 个点对应的次数最小的多项式。这里假设 \(x_i, y_i\) 分别互不相同。如果能对每个点 \(i\) 找到一个函数 \(f_i(x)\)，使得 \(f_i(x_i) = 1\)，并且 \(\forall j\neq i,\ f_i(x_j) = 0\)，那我们直接对所有点值求和得到

\[F(x) = \sum_{i=1}^n y_i \times f_{i}(x) \]

这样 \(F(x)\) 肯定经过这 \(n\) 个点。随后我们就需要满足一下次数的最小性了。可以发现我们需要的 \(f_i(x)\) 应当是 \(n-1\) 次多项式，因此这启发我们通过 \(n-1\) 个单项式相乘的方法去限制 \(x = x_j\) 时 \(f_i(x)\) 为 \(0\)，直接乘入一个单项式 \((x - x_j)\) 即可，这就得到了 \(f_i(x)\) 的部分构造。然后我们发现还差一个 \(f_i(x_i) = 1\)，这如何满足呢？考虑一个 \(p / p = 1\)，我们将目前构造出的多项式除以 \(\prod (x_i - x_j)\)，这样在 \(x = x_i\) 时上下相同，取到 \(1\)。
很容易验证这样得出的 \(f_i(x)\) 是满足所有要求的。因此我们构造出了 \(F(x)\) 的一般形式。

总结一下能得到拉格朗日插值法的公式：

\[F(x) = \sum_{i=1}^n y_i \prod_{j\neq i}\frac{x - x_j}{x_i - x_j} \]

实现

inline int Lagrange(int n, int x[], int y[], int x_qry) {
	int ret = 0;
	rep(i,0,n) {
		int s1 = 1, s2 = 1;
		rep(j,0,n) {
			if (i != j) {
				s1 = 1ll * s1 * (x_qry - x[j]) % mod;
				s2 = 1ll * s2 * (x[i] - x[j]) % mod;
			} 
		} ret = (ret + 1ll * y[i] * s1 % mod * qp(s2, mod - 2) % mod) % mod;
	} return (ret + mod) % mod; 
} 

由这个式子可以得到一种 \(x_i\) 点值连续时的 \(O(n)\) 做法。不展开，留作练习。

实现

int s1[N], s2[N], inv[N];
inline int CLagrange(int n, int x[], int y[], int x_qry) {
	int ret = 0;
	s1[0] = (x_qry - x[0]) % mod, s2[n+1] = 1;
	rep(i,1,n) s1[i] = 1ll * s1[i-1] * (x_qry - x[i]) % mod;
	pre(i,n,0) s2[i] = 1ll * s2[i+1] * (x_qry - x[i]) % mod;
	inv[0] = inv[1] = 1;
	rep(i,2,n) inv[i] = -1ll * mod / i * inv[mod % i] % mod;
	rep(i,2,n) inv[i] = 1ll * inv[i-1] * inv[i] % mod;
	rep(i,0,n) ret = (ret + (1ll * y[i] * (i == 0 ? 1 : s1[i-1]) % mod * s2[i+1] % mod * inv[i] % mod * (((n-i) & 1) ? -1 : 1) * inv[n-i]) % mod) % mod;
	return (ret + mod) % mod;
} 

然后有一种动态做这个东西的方法，叫重心拉格朗日插值法，加入一个点的复杂度是 \(O(n)\) 的。

具体的式子不再推了，就是提取出一个 \(\prod_{i} (x - x_i)\)，最后的部分除了形如 \((x_i - x_j)\) 的系数就只剩下 \(y_i / (x - x_i)\) 了。我们发现把系数组合得到 \(t_i = y_i / \prod_{j\neq i} (x_i - x_j)\) 后就能得到一个每个点计算贡献为 \(O(n)\) 的式子。

我们容易实现出 \(O(n^2)\) 得到所有系数的程序。如下：

实现

int qp(int a, int b) {
    int ret = 1;
    while (b) {
        if (b & 1) ret = 1ll * ret * a % mod;
        a = 1ll * a * a % mod;
        b >>= 1;
    } return ret;
}

struct poly : vector<int> {
    int operator() (const int& x) {
        int ret = 0;
        for (int i = 0, p = 1; i < size(); ++ i, p = 1ll * p * x % mod)  {
            ret = (ret + 1ll * p * (*this)[i]) % mod;
        } return ret;
    }
};
int n, k, x[N], y[N], inv[N];
poly f, sum, tmp;
void lag(int n, int x[], int y[]) { 
    sum.resize(n + 1), tmp.resize(n + 1); sum[0] = 1;
    for (int i = 1; i <= n; ++ i, swap(sum, tmp)) {
        tmp[0] = 0, inv[i] = qp(mod - x[i], mod - 2);
        rep(j,1,i) tmp[j] = sum[j - 1];
        rep(j,0,i) tmp[j] = (tmp[j] + 1ll * sum[j] * (mod - x[i])) % mod;
    } f.resize(n + 1);
    for (int i = 1, iv; i <= n; ++ i) {
        iv = 1;
        rep(j,1,n) if (j != i) 
            iv = 1ll * iv * (x[i] - x[j] + mod) % mod;
        iv = 1ll * qp(iv, mod - 2) * y[i] % mod;
        for (int j = 0, lst = 0; j < n; ++ j) {
            tmp[j] = 1ll * (sum[j] - lst + mod) * inv[i] % mod;
            f[j] = (f[j] + 1ll * iv * tmp[j]) % mod, lst = tmp[j];
        }
    } 
}

一个具体的应用是求

\[S_k(n) = \sum_{i=1}^n i^k \]

我们可以通过作差的方法得到 \(S_k(n)\) 是一个关于 \(n\) 的 \(k + 1\) 阶多项式。
对 \(S_k(n)\) 做差得到了 \(n^k\)，这个东西是关于 \(n\) 的 \(k\) 阶多项式，因此 \(S_k(n)\) 是一个关于 \(n\) 的 \(k + 1\) 阶多项式。

随后可以直接通过拉格朗日差值得到系数和点值。

我们通过上面的讨论得知在点值连续时可以 \(O(k)\) 地得到这个多项式，这里只需要计算点值所以容易些。对于 \(i\in [1, k + 2]\) 计算 \(i^{k}\) 是容易 \(O(k)\) 的，于是我们就得到了一种 \(O(k)\) 计算 \(S_k(n)\) 的方法。

实现

int S(int n, int k) {
    vector<int> pref(k + 4), suff(k + 4), pw(k + 4), prime(k + 4); 
    vector<bool> __vis(k + 4, 0);
    pw[1] = 1;
    for (int i = 2, cnt = 0; i <= k + 2; ++ i) {
        if (!__vis[i]) prime[++ cnt] = i, pw[i] = qp(i, k);
        for (int j = 1; j <= cnt and i * prime[j] <= k + 2; ++ j) {
            __vis[i * prime[j]] = 1;
            pw[i * prime[j]] = mul(pw[i], pw[prime[j]]);
            if (i % prime[j] == 0) break;
        }
    } rep(i,2,k+2) pw[i] = add(pw[i], pw[i - 1]);
    pref[0] = suff[k + 3] = 1;
    rep(i, 1, k + 2) pref[i] = mul(pref[i - 1], (n - i + mod));
    pre(i, k + 2, 1) suff[i] = mul(suff[i + 1], (n - i + mod));
    int tmp = 0, ret = 0;
    rep(i, 1, k + 2)
        ret = add(ret, mul(((k - i) & 1) ? mod - 1 : 1, pw[i], pref[i - 1], suff[i + 1], ifc(i - 1), ifc(k + 2 - i)));
    return ret;
}

同样的思路似乎可以自然地应用到多元/高维拉格朗日插值。