我在做 topcoder13369,大体思路是:矩阵每个元素放一个一次多项式,最后矩阵行列式是一个 n−1 次多项式,给 x 代入 n 个值进去高斯消元得到 n 个点值,再拉格朗日插值得到多项式系数。
问题出在,一开始我 x 从 0 开始代入值,一直不对,改成 x 从 1 开始代入值就对了。
不太懂为啥拉格朗日插值为什么不能处理 x=0 的点。也可能是别的地方的锅?
topcoder13369 的代码
#include <bits/stdc++.h>
using i64 = long long;
const int P = 1e9 + 7;
void inc(int &x, const int y) {
x += y;
if (x >= P) x -= P;
}
int plus(const int x, const int y) {
return (x + y) >= P ? (x + y - P) : (x + y);
}
void mul(int &x, const int y) {
x = (i64) x * y % P;
}
int times(const int x, const int y) {
return (i64) x * y % P;
}
int power(int a, int b) {
int res = 1;
while (b) {
if (b & 1) mul(res, a);
mul(a, a);
b >>= 1;
}
return res;
}
int inverse(int a) {
return power(a, P - 2);
}
class TreeDistance {
private:
int guass(int n, std::vector<std::vector<int>> &a) {
int det = 1;
for (int i = 0; i < n; ++i) {
if (!a[i][i]) {
for (int j = i + 1; j < n; ++j) {
if (a[j][i]) {
std::swap(a[i], a[j]);
mul(det, P - 1);
break;
}
}
}
if (!a[i][i])
return 0;
mul(det, a[i][i]);
int inv = inverse(a[i][i]);
for (int j = i; j < n; ++j)
mul(a[i][j], inv);
for (int j = i + 1; j < n; ++j) {
int d = a[j][i];
for (int k = i; k < n; ++k)
inc(a[j][k], P - times(a[i][k], d));
}
}
return det;
}
std::vector<int> lagrange(int n, std::vector<int> x, std::vector<int> y) {
std::vector<int> l(n + 1), t(n + 1), res(n);
l[0] = 1;
for (int i = 0; i < n; ++i) {
for (int j = 1; j <= n; ++j)
t[j] = l[j - 1];
t[0] = 0;
for (int j = 0; j <= n; ++j)
inc(t[j], times(l[j], P - x[i]));
std::swap(l, t);
}
for (int i = 0; i < n; ++i) {
int p = 1;
for (int j = 0; j < n; ++j)
if (j != i)
mul(p, plus(x[i], P - x[j]));
p = times(y[i], inverse(p));
int inv = inverse(P - x[i]);
int last = 0;
for (int j = 0; j < n; ++j) {
int now = times(plus(l[j], P - last), inv);
inc(res[j], times(now, p));
last = now;
}
}
return res;
}
public:
int countTrees(std::vector<int> p, int k) {
int n = p.size() + 1;
std::vector<std::vector<std::pair<int, int>>> mat(n, std::vector<std::pair<int, int>>(n, {P - 1, 0}));
for (int i = 0; i < n; ++i)
mat[i][i].first = n - 1;
for (int i = 1; i < n; ++i) {
int j = p[i - 1];
inc(mat[i][j].first, 1), inc(mat[i][j].second, P - 1);
inc(mat[j][i].first, 1), inc(mat[j][i].second, P - 1);
inc(mat[i][i].first, P - 1), inc(mat[i][i].second, 1);
inc(mat[j][j].first, P - 1), inc(mat[j][j].second, 1);
}
// 代值过程在这里
std::vector<int> fx(n), fy(n);
for (int x = 0; x < n; ++x) {
std::vector<std::vector<int>> a(n - 1, std::vector<int>(n - 1));
for (int i = 0; i < n - 1; ++i)
for (int j = 0; j < n - 1; ++j)
a[i][j] = (mat[i][j].first + times(mat[i][j].second, x)) % P;
fx[x] = x;
fy[x] = guass(n - 1, a);
}
auto f = lagrange(n, fx, fy);
int ans = 0;
for (int i = std::max(n - k - 1, 0); i < n; ++i)
inc(ans, f[i]);
return ans;
}
};
int main() {
std::cout << TreeDistance().countTrees(
{0, 1, 2, 2, 0}, 2
) << '\n';
return 0;
}
真的是萌新,拉格朗日插值就背了个式子,不太懂原理。