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)$.

Hình trên là đa thức đi qua các điểm đã cho; hình dưới là đa thức cơ sở $L_i$, bằng 1 tại nút $i$ và 0 tại các nút khác. Đổi $y$ hoặc $i$ để xem.

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;
#define 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';
}

ChạyInput2 4
1 5
2 7

Output11

Độ phức tạp: $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

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;
#define 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';
}

ChạyInput2
1 4 9
3

Output16

Độ phức tạp: $O(N + \log p)$.

Tổng lũy thừa

$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 nút liên tiếp. 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;
#define 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';
}

ChạyInput4 1

Output10

$O(k \log k)$ cho các lũy thừa, $O(k)$ cho nội suy. Phần lũy thừa là chỗ tốn nhất, và làm nhanh hơn được vì $i^k$ là hàm nhân tính hoàn toàn: chỉ cần lũy thừa nhanh tại các số nguyên tố, rồi suy ra các số khác bằng một phép nhân qua sàng tuyến tính.

Sàng tuyến tính nhanh hơn lũy thừa nhanh từng số khoảng 3 lần. Quy tắc: cần $i^k$ cho mọi $i$ liên tiếp thì chỉ tính tại số nguyên tố.

DP 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 mọi DP có 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;
#define 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';
}

ChạyInput1 4

Output4

Độ phức tạp: quy hoạch động $O(N^2)$, nội suy $O(N)$.