Lagrange Interpolation

Đa thức cơ sở

Với các nút phân biệt $x_1, \dots, x_n$ và giá trị $y_i$, đặt

$L_i$ có bậc $n - 1$, bằng 1 tại $x_i$ và bằng 0 tại mọi nút khác, nên $f(x) = \sum_i y_i L_i(x)$ đi qua đủ $n$ điểm. Hai đa thức bậc nhỏ hơn $n$ trùng nhau tại $n$ điểm thì hiệu của chúng có $n$ nghiệm nên bằng 0: đa thức nội suy là duy nhất.

Theo mô-đun nguyên tố $p$, mọi phép chia thành nhân nghịch đảo, chỉ cần các nút phân biệt theo mô-đun $p$. Với mỗi $i$, nhân tử số và mẫu số riêng rồi lấy một nghịch đảo, tổng $O(n^2 + n \log p)$.

Ví dụGiá Trị Qua Các Điểm

Đề. $n \le 2000$ điểm $(x_i, y_i)$ với $x_i$ phân biệt, tính $f(k)$ theo mô-đun $998244353$ với $f$ là đa thức bậc nhỏ hơn $n$ qua các điểm đó.

Lời giải. Áp thẳng công thức nội suy tại điểm $k$:

Với mỗi $i$, nhân tử số và mẫu số riêng trong một vòng $O(n)$, rồi lấy một nghịch đảo của mẫu số.

num bắt đầu bằng y[i] để nhân thẳng vào tử số; den là tích các hiệu $x_i - x_j$. Hiệu có thể âm, mà phép % của C++ giữ dấu, nên đưa về $[0, p)$ trước khi nhân. Các $x_i$ phân biệt nên mẫu số khác 0 và có nghịch đảo.

for (int i = 0; i < n; i++) {
    ll num = y[i], den = 1;
    for (int j = 0; j < n; j++)
        if (j != i) {
            num = num * ((k - x[j]) % MOD + MOD) % MOD;
            den = den * ((x[i] - x[j]) % MOD + MOD) % MOD;
        }
    ans = (ans + num * power(den, MOD - 2)) % MOD;
}

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const ll MOD = 998244353;
ll n, k, x[2005], y[2005];

ll power(ll a, ll b) {
    ll r = 1;
    a %= MOD;
    while (b) {
        if (b & 1) r = r * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    return r;
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n >> k;
    for (int i = 0; i < n; i++) cin >> x[i] >> y[i];
    ll ans = 0;
    for (int i = 0; i < n; i++) {
        ll num = y[i], den = 1;
        for (int j = 0; j < n; j++)
            if (j != i) {
                num = num * ((k - x[j]) % MOD + MOD) % MOD;
                den = den * ((x[i] - x[j]) % MOD + MOD) % MOD;
            }
        ans = (ans + num * power(den, MOD - 2)) % MOD;
    }
    cout << ans << '\n';
}

Input2 4
1 5
2 7

Output11

$O(n^2 + n \log p)$, khoảng $4 \cdot 10^6$ phép nhân với $n = 2000$.

Nút liên tiếp trong O(n)

Khi các nút là $0, 1, \dots, n$, mẫu số có dạng đóng: $\prod_{j \ne i} (i - j) = i! \cdot (n - i)! \cdot (-1)^{n - i}$. Tử số $\prod_{j \ne i}(T - j)$ là tích tiền tố $\prod_{j < i}(T - j)$ nhân tích hậu tố $\prod_{j > i}(T - j)$. Tính sẵn giai thừa và nghịch đảo giai thừa, mỗi số hạng $O(1)$, tổng $O(n)$.

Nếu $T \le n$ thì trả ngay $P(T)$, vì khi đó tử số của mọi cơ sở trừ một đều bằng 0.

Ví dụBảng Giá Trị Liên Tiếp

Đề. Đa thức bậc không quá $N \le 10^5$, biết $P(0), \dots, P(N)$ theo mô-đun $10^9 + 7$. Tính $P(T)$ với $T \le 10^9$.

Lời giải. Nút là $0, 1, \dots, N$ liên tiếp. Mẫu số của cơ sở $i$ là $\prod_{j \ne i}(i - j) = i! \cdot (N - i)! \cdot (-1)^{N - i}$: các $j < i$ cho tích $i!$, các $j > i$ cho $(N - i)!$ với $N - i$ dấu trừ. Tử số $\prod_{j \ne i}(T - j)$ là tích tiền tố $\prod_{j < i}$ nhân tích hậu tố $\prod_{j > i}$. Tính sẵn hai mảng tích và bảng nghịch đảo giai thừa thì mỗi số hạng $O(1)$.

Bước 1: trường hợp có sẵn. $T \le N$ thì trả ngay a[t]. Trường hợp này cũng tránh việc tử số có một thừa số bằng 0.

Bước 2: giai thừa nghịch đảo. Một lần lũy thừa cho ifac[n], rồi lùi dần: $\frac{1}{(i-1)!} = \frac{i}{i!}$.

fac[0] = 1;
for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i % MOD;
ifac[n] = power(fac[n], MOD - 2);
for (int i = n; i >= 1; i--) ifac[i - 1] = ifac[i] * i % MOD;

Bước 3: tích tiền tố và hậu tố. pre[i] $= \prod_{j < i}(T - j)$, suf[i] $= \prod_{j \ge i}(T - j)$; tử số của cơ sở $i$ là pre[i] * suf[i + 1].

pre[0] = 1;
for (int i = 1; i <= n + 1; i++) pre[i] = pre[i - 1] * ((t - (i - 1)) % MOD) % MOD;
suf[n + 1] = 1;
for (int i = n; i >= 0; i--) suf[i] = suf[i + 1] * ((t - i) % MOD) % MOD;

Bước 4: cộng các số hạng. Dấu $(-1)^{N - i}$ xử lý bằng cách cộng hay trừ.

for (int i = 0; i <= n; i++) {
    ll term = a[i] * pre[i] % MOD * suf[i + 1] % MOD * ifac[i] % MOD * ifac[n - i] % MOD;
    if ((n - i) % 2) ans = (ans - term + MOD) % MOD;
    else ans = (ans + term) % MOD;
}

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const ll MOD = 1e9 + 7;
const int N = 1e5 + 5;
ll n, t, a[N], pre[N], suf[N], fac[N], ifac[N];

ll power(ll x, ll b) {
    ll r = 1;
    x %= MOD;
    while (b) {
        if (b & 1) r = r * x % MOD;
        x = x * x % MOD;
        b >>= 1;
    }
    return r;
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n;
    for (int i = 0; i <= n; i++) cin >> a[i];
    cin >> t;
    if (t <= n) {
        cout << a[t] << '\n';
        return 0;
    }
    fac[0] = 1;
    for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i % MOD;
    ifac[n] = power(fac[n], MOD - 2);
    for (int i = n; i >= 1; i--) ifac[i - 1] = ifac[i] * i % MOD;
    pre[0] = 1;
    for (int i = 1; i <= n + 1; i++) pre[i] = pre[i - 1] * ((t - (i - 1)) % MOD) % MOD;
    suf[n + 1] = 1;
    for (int i = n; i >= 0; i--) suf[i] = suf[i + 1] * ((t - i) % MOD) % MOD;
    ll ans = 0;
    for (int i = 0; i <= n; i++) {
        ll term = a[i] * pre[i] % MOD * suf[i + 1] % MOD * ifac[i] % MOD * ifac[n - i] % MOD;
        if ((n - i) % 2) ans = (ans - term + MOD) % MOD;
        else ans = (ans + term) % MOD;
    }
    cout << ans << '\n';
}

Input2
1 4 9
3

Output16

$O(N + \log p)$.

Tổng lũy thừa là đa thức

$S(n) = \sum_{i=1}^{n} i^k$ là đa thức bậc $k + 1$ theo $n$: hiệu $S(n) - S(n - 1) = n^k$ có bậc $k$, và một dãy có hiệu bậc $k$ thì chính nó có bậc $k + 1$ (như tích phân tăng bậc). Vậy chỉ cần $S(0), S(1), \dots, S(k + 1)$, rồi nội suy tại $n$ bằng cách nút liên tiếp, tổng $O(k \log k)$ cho việc tính các lũy thừa và $O(k)$ cho nội suy. Không cần công thức Faulhaber hay số Bernoulli.

Cách nhận ra “là đa thức” nói chung: lấy hiệu liên tiếp nhiều lần; nếu đáp án bậc $d$ thì hiệu lần thứ $d + 1$ bằng 0. Viết brute force in thử vài giá trị là kiểm được nhanh.

Ví dụTổng Các Lũy Thừa

Đề. Tính $1^k + 2^k + \dots + n^k$ theo mô-đun $10^9 + 7$ với $n \le 10^9$, $k \le 10^6$.

Lời giải. $S(n) = \sum_{i=1}^n i^k$ là đa thức bậc $k + 1$ theo $n$, nên xác định bởi $k + 2$ giá trị $S(0), \dots, S(k + 1)$. Tính các giá trị đó bằng cộng dồn $i^k$, rồi nội suy tại $n$ bằng đúng khung nút liên tiếp của ví dụ trước với $m = k + 1$ thay cho $N$.

Bước 1: các giá trị tại nút. s[i] cộng dồn power(i, k); đây là phần chậm nhất. n <= m thì s[n] đã là đáp án.

int m = k + 1;
for (int i = 1; i <= m; i++) s[i] = (s[i - 1] + power(i, k)) % MOD;
if (n <= m) {
    cout << s[n] << '\n';
    return 0;
}

Bước 2: nội suy chép nguyên khung của ví dụ trước: pre, suf, giai thừa nghịch đảo và dấu $(-1)^{m - i}$.

Kiểm nhanh với mẫu: $k = 1$, $n = 4$ ra $1 + 2 + 3 + 4 = 10$.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const ll MOD = 1e9 + 7;
const int N = 1e6 + 5;
ll n, k, s[N], pre[N], suf[N], fac[N], ifac[N];

ll power(ll x, ll b) {
    ll r = 1;
    x %= MOD;
    while (b) {
        if (b & 1) r = r * x % MOD;
        x = x * x % MOD;
        b >>= 1;
    }
    return r;
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n >> k;
    int m = k + 1;
    for (int i = 1; i <= m; i++) s[i] = (s[i - 1] + power(i, k)) % MOD;
    if (n <= m) {
        cout << s[n] << '\n';
        return 0;
    }
    fac[0] = 1;
    for (int i = 1; i <= m; i++) fac[i] = fac[i - 1] * i % MOD;
    ifac[m] = power(fac[m], MOD - 2);
    for (int i = m; i >= 1; i--) ifac[i - 1] = ifac[i] * i % MOD;
    pre[0] = 1;
    for (int i = 1; i <= m + 1; i++) pre[i] = pre[i - 1] * ((n - (i - 1)) % MOD) % MOD;
    suf[m + 1] = 1;
    for (int i = m; i >= 0; i--) suf[i] = suf[i + 1] * ((n - i) % MOD) % MOD;
    ll ans = 0;
    for (int i = 0; i <= m; i++) {
        ll term = s[i] * pre[i] % MOD * suf[i + 1] % MOD * ifac[i] % MOD * ifac[m - i] % MOD;
        if ((m - i) % 2) ans = (ans - term + MOD) % MOD;
        else ans = (ans + term) % MOD;
    }
    cout << ans << '\n';
}

Input4 1

Output10

$O(k \log k)$ cho các lũy thừa, $O(k)$ cho nội suy.

DP có tham số lớn

Có những DP mà kết quả, xét như hàm của một tham số lớn $D$, là đa thức bậc nhỏ. Ví dụ đếm cách gán nhãn trong $[1, D]$ cho cây sao cho con không lớn hơn cha: gọi $G_v(x)$ là số cách gán cây con $v$ với nhãn của $v$ không vượt $x$. Lá có $G_v(x) = x$, bậc 1. Ở đỉnh trong, số cách khi nhãn của $v$ đúng bằng $y$ là $\prod_c G_c(y)$, có bậc $sz_v - 1$; cộng dồn theo $y \le x$ tăng bậc lên 1. Vậy $G_v$ có bậc $sz_v$, và $G_{\text{gốc}}$ có bậc $n$.

Tính $G$ cho $x = 0..n$ bằng DP $O(n^2)$, rồi nội suy $G_{\text{gốc}}(D)$ với $D$ tới $10^9$. Cùng mẹo cho DP mà một chiều là “giá trị tối đa $D$”: chạy DP với $D$ nhỏ, chứng minh bậc, nội suy.

Ví dụNhãn Không Tăng Trên Cây

Đề. Cây $N \le 3000$ đỉnh, cha của $v$ là $p_v < v$. Đếm cách gán nhãn trong $[1, D]$, $D \le 10^9$, sao cho nhãn con không lớn hơn nhãn cha, theo mô-đun $10^9 + 7$.

Lời giải. Như phân tích ở trên, $G_v(x)$ là số cách gán cây con của $v$ khi nhãn của $v$ không vượt $x$, có bậc $sz_v$ theo $x$; gốc cho đa thức bậc $N$. Tính $G$ tại $x = 0..N$ bằng quy hoạch động rồi nội suy tại $D$ bằng khung nút liên tiếp. Vì $p_v < v$, duyệt $v$ từ $N$ về 1 thì mọi con được xử lý trước cha.

Bước 1: nhân các con. f[v][x] ban đầu là 1 với mọi $x$. Khi con $c$ xong, f[par[c]][x] nhân thêm $G_c(x)$. Tới lượt $v$, mọi con đã xong, nên f[v][y] $= \prod_c G_c(y)$: số cách khi nhãn của $v$ đúng bằng $y$.

Bước 2: đổi thành tổng tiền tố. Đặt f[v][0] = 0 (nhãn phải từ 1), rồi cộng dồn: giờ f[v][x] $= \sum_{y \le x} \prod_c G_c(y) = G_v(x)$.

for (int v = n; v >= 1; v--) {
    f[v][0] = 0;
    for (int x = 1; x <= n; x++) f[v][x] = (f[v][x - 1] + f[v][x]) % MOD;
    if (v > 1)
        for (int x = 0; x <= n; x++) f[par[v]][x] = (ll)f[par[v]][x] * f[v][x] % MOD;
}

Bước 3: nội suy. d <= n thì f[1][d] đã là đáp án; ngược lại nội suy f[1][0..n] tại $d$. Bảng f là $3001 \times 3001$ số int, khoảng 36 MB.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

const ll MOD = 1e9 + 7;
const int N = 3005;
int n, par[N], f[N][N];
ll d, pre[N], suf[N], fac[N], ifac[N];

ll power(ll x, ll b) {
    ll r = 1;
    x %= MOD;
    while (b) {
        if (b & 1) r = r * x % MOD;
        x = x * x % MOD;
        b >>= 1;
    }
    return r;
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n >> d;
    for (int v = 2; v <= n; v++) cin >> par[v];
    for (int v = 1; v <= n; v++)
        for (int x = 0; x <= n; x++) f[v][x] = 1;
    for (int v = n; v >= 1; v--) {
        f[v][0] = 0;
        for (int x = 1; x <= n; x++) f[v][x] = (f[v][x - 1] + f[v][x]) % MOD;
        if (v > 1)
            for (int x = 0; x <= n; x++) f[par[v]][x] = (ll)f[par[v]][x] * f[v][x] % MOD;
    }
    if (d <= n) {
        cout << f[1][d] << '\n';
        return 0;
    }
    fac[0] = 1;
    for (int i = 1; i <= n; i++) fac[i] = fac[i - 1] * i % MOD;
    ifac[n] = power(fac[n], MOD - 2);
    for (int i = n; i >= 1; i--) ifac[i - 1] = ifac[i] * i % MOD;
    pre[0] = 1;
    for (int i = 1; i <= n + 1; i++) pre[i] = pre[i - 1] * ((d - (i - 1)) % MOD) % MOD;
    suf[n + 1] = 1;
    for (int i = n; i >= 0; i--) suf[i] = suf[i + 1] * ((d - i) % MOD) % MOD;
    ll ans = 0;
    for (int i = 0; i <= n; i++) {
        ll term = f[1][i] * pre[i] % MOD * suf[i + 1] % MOD * ifac[i] % MOD * ifac[n - i] % MOD;
        if ((n - i) % 2) ans = (ans - term + MOD) % MOD;
        else ans = (ans + term) % MOD;
    }
    cout << ans << '\n';
}

Input1 4

Output4

Quy hoạch động $O(N^2)$, nội suy $O(N)$.