跳转至

拉格朗日插值

\(f(x)=\sum_{i=1}^ny_i*\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}\)

\(k\) 阶多项式需要 \(k+1\) 个点插

例题

P4781 【模板】拉格朗日插值 - 洛谷

代码
C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
int __OneWan_2024 = [](){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    return 0;
}();
const i64 P = 998244353;
i64 qpow(i64 a, i64 b, i64 res = 1) {
    while (b) {
        if (b & 1) {
            res = res * a % P;
        }
        a = a * a % P;
        b >>= 1;
    }
    return res;
}
i64 x[2005], y[2005];
int main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    int n, k;
    cin >> n >> k;
    for (int i = 0 ; i < n ; i++) {
        cin >> x[i] >> y[i];
    }
    auto larange = [&](long long k) {
        long long res = 0;
        for (int i = 0 ; i < n ; i++) {
            i64 s1 = 1, s2 = 1;
            for (int j = 0 ; j < n ; j++) {
                if (i == j) continue;
                s1 = s1 * (k - x[j]) % P;
                s2 = s2 * (x[i] - x[j]) % P;
            }
            res = (res + y[i] * s1 % P * qpow(s2, P - 2) % P + P) % P;
        }
        return res;
    };
    cout << larange(k) << "\n";
    return 0;
}
习题

P4593 [TJOI2018] 教科书般的亵渎 - 洛谷

\(ans=\sum_{i=0}^m\left(\sum_{j=1}^{n-a_i}j^{m+1}-\sum_{j=i+1}^m({(a_j-a_i)}^{m+1})\right)\)

代码
C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
int __OneWan_2024 = [](){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    return 0;
}();
const i64 P = 1000000007;
i64 qpow(i64 a, i64 b, i64 res = 1) {
    while (b) {
        if (b & 1) {
            res = res * a % P;
        }
        a = a * a % P;
        b >>= 1;
    }
    return res;
}
i64 y[55];
i64 a[55];
void solve() {
    int n, m;
    cin >> n >> m;
    for (int i = 1 ; i <= m ; i++) {
        cin >> a[i];
    }
    sort(a, a + m + 1);
    int k = m + 1;
    y[1] = 1;
    for (int i = 2 ; i <= k + 2 ; i++) {
        y[i] = (y[i - 1] + qpow(i, k)) % P;
    }
    auto larange = [&](i64 x) {
        if (x <= k + 2) return y[x];
        i64 res = 0;
        for (int i = 1 ; i <= k + 2 ; i++) {
            i64 s1 = 1, s2 = 1;
            for (int j = 1 ; j <= k + 2 ; j++) {
                if (i == j) continue;
                s1 = s1 * (x - j) % P;
                s2 = s2 * (i - j) % P;
            }
            res = (res + y[i] * s1 % P * qpow(s2, P - 2) % P + P) % P;
        }
        return res;
    };
    i64 ans = 0;
    for (int i = 0 ; i <= m ; i++) {
        ans = (ans + larange(n - a[i])) % P;
        i64 res = 0;
        for (int j = i + 1 ; j <= m ; j++) {
            res = res + qpow(a[j] - a[i], k);
        }
        ans = ((ans - res) % P + P) % P;
    }
    cout << ans << "\n";
}
int main() {
    int T;
    cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}

Problem - F - Codeforces

求 \(\sum_{i=1}^ni^k\)

代码
C++
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
int __OneWan_2024 = [](){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    return 0;
}();
const i64 P = 1000000007;
i64 pre[1000005], suf[1000005], y[1000005], fact[1000005], invf[1000005];
i64 qpow(i64 a, i64 b = P - 2, i64 res = 1) {
    while (b) {
        if (b & 1) {
            res = res * a % P;
        }
        a = a * a % P;
        b >>= 1;
    }
    return res;
}
int main() {
    int n, k;
    cin >> n >> k;
    for (int i = 1 ; i <= k + 2 ; i++) {
        y[i] = (y[i - 1] + qpow(i, k)) % P;
    }
    fact[0] = 1;
    for (int i = 1 ; i <= k + 2 ; i++) {
        fact[i] = fact[i - 1] * i % P;
    }
    invf[k + 2] = qpow(fact[k + 2], P - 2);
    for (int i = k + 2 ; i >= 1 ; i--) {
        invf[i - 1] = invf[i] * i % P;
    }
    auto lagrange = [&](int x) {
        if (x <= k + 2) {
            return y[x];
        }
        pre[0] = 1;
        for (int i = 1 ; i <= k + 2 ; i++) {
            pre[i] = pre[i - 1] * (x - i) % P;
        }
        suf[k + 3] = 1;
        for (int i = k + 2 ; i >= 1 ; i--) {
            suf[i] = suf[i + 1] * (x - i) % P;
        }
        // 分子 pre[i - 1] * suf[i + 1]
        // 分母 invf[i - 1] * invf[k + 2 - i] 符号(-1)^{k + 2 - i}
        i64 res = 0;
        for (int i = 1 ; i <= k + 2 ; i++) {
            i64 s1 = pre[i - 1] * suf[i + 1] % P;
            i64 s2 = invf[i - 1] * invf[k + 2 - i] % P;
            if ((k + 2 - i) & 1) {
                s2 = P - s2;
            }
            res = (res + s1 * s2 % P * y[i] % P) % P;
        }
        return res;
    };
    cout << lagrange(n);
    return 0;
}

tyvj 1858 XLkxc - HydroOJ

f函数是k+1阶多项式, g函数是k+2阶多项式

设答案为h函数, 则 h函数是k+3阶多项式

代码
C++
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
#include <bits/stdc++.h>
using namespace std;
using i64 = long long;
int __OneWan_2024 = [](){
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    return 0;
}();
const i64 P = 1234567891;
i64 pre[200], suf[200], y[200], fact[200], invf[200];
i64 qpow(i64 a, i64 b = P - 2) {
    i64 res = 1;
    while (b) {
        if (b & 1) {
            res = res * a % P;
        }
        a = a * a % P;
        b >>= 1;
    }
    return res;
}
i64 f[200], g[200], h[200];
void solve() {
    int k, a, n, d;
    cin >> k >> a >> n >> d;
    for (int i = 1 ; i <= k + 3 ; i++) {
        f[i] = (f[i - 1] + qpow(i, k)) % P;
        g[i] = (g[i - 1] + f[i]) % P;
    }
    auto lagrange1 = [&](i64 x) {
        if (x <= k + 3) {
            return g[x];
        }
        pre[0] = 1;
        for (int i = 1 ; i <= k + 3 ; i++) {
            pre[i] = pre[i - 1] * ((x - i) % P) % P;
        }
        suf[k + 4] = 1;
        for (int i = k + 3 ; i >= 1 ; i--) {
            suf[i] = suf[i + 1] * ((x - i) % P) % P;
        }
        // 分子 pre[i - 1] * suf[i + 1]
        // 分母 invf[i - 1] * invf[k + 3 - i] 符号(-1)^{k + 3 - i}
        i64 res = 0;
        for (int i = 1 ; i <= k + 3 ; i++) {
            i64 s1 = pre[i - 1] * suf[i + 1] % P;
            i64 s2 = invf[i - 1] * invf[k + 3 - i] % P;
            if ((k + 3 - i) & 1) {
                s2 = P - s2;
            }
            res = (res + s1 * s2 % P * g[i] % P) % P;
        }
        return res;
    };
    auto lagrange2 = [&](i64 x) {
        if (x <= k + 4) {
            return h[x];
        }
        pre[0] = 1;
        for (int i = 1 ; i <= k + 4 ; i++) {
            pre[i] = pre[i - 1] * ((x - i) % P) % P;
        }
        suf[k + 5] = 1;
        for (int i = k + 4 ; i >= 1 ; i--) {
            suf[i] = suf[i + 1] * ((x - i) % P) % P;
        }
        // 分子 pre[i - 1] * suf[i + 1]
        // 分母 invf[i - 1] * invf[k + 4 - i] 符号(-1)^{k + 4 - i}
        i64 res = 0;
        for (int i = 1 ; i <= k + 4 ; i++) {
            i64 s1 = pre[i - 1] * suf[i + 1] % P;
            i64 s2 = invf[i - 1] * invf[k + 4 - i] % P;
            if ((k + 4 - i) & 1) {
                s2 = P - s2;
            }
            res = (res + s1 * s2 % P * h[i] % P) % P;
        }
        return res;
    };
    for (int i = 1 ; i <= k + 4 ; i++) {
        h[i] = (h[i - 1] + lagrange1(a + 1LL * (i - 1) * d)) % P;
    }
    cout << lagrange2(n + 1) << "\n";
}

int main() {
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    fact[0] = 1;
    for (int i = 1 ; i <= 150 ; i++) {
        fact[i] = fact[i - 1] * i % P;
    }
    invf[150] = qpow(fact[150], P - 2);
    for (int i = 150 ; i >= 1 ; i--) {
        invf[i - 1] = invf[i] * i % P;
    }
    int T;
    cin >> T;
    while (T--) {
        solve();
    }
    return 0;
}