拉格朗日插值为什么不能处理 x=0 的点?
  • 板块学术版
  • 楼主RenaMoe
  • 当前回复1
  • 已保存回复1
  • 发布时间2021/7/23 09:17
  • 上次更新2023/11/4 13:36:34
查看原帖
拉格朗日插值为什么不能处理 x=0 的点?
52243
RenaMoe楼主2021/7/23 09:17

我在做 topcoder13369,大体思路是:矩阵每个元素放一个一次多项式,最后矩阵行列式是一个 n1n-1 次多项式,给 xx 代入 nn 个值进去高斯消元得到 nn 个点值,再拉格朗日插值得到多项式系数。

问题出在,一开始我 xx00 开始代入值,一直不对,改成 xx11 开始代入值就对了。

不太懂为啥拉格朗日插值为什么不能处理 x=0x=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;
}

真的是萌新,拉格朗日插值就背了个式子,不太懂原理。

2021/7/23 09:17
加载中...