闲话
数学这里证明是个避不过去的东西。你总得面对他。
—— 巨佬
今日 ABC(
先写一点 慢慢补上
拉格朗日插值
熟知平面上 \(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)\) 的式子。
一个具体的应用是求
\[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;
}
未完。
标签:int,闲话,inv,ret,23.1,s2,s1,mod From: https://www.cnblogs.com/joke3579/p/chitchat230107.html