Trong một đề luyện tập Tin học trẻ của trường THPT Kim Liên vừa rồi có bài toán sau:

Cho ~n, k~ (~1 \leq n \leq 10^8, 1 \leq k \leq 10^3~), tính giá trị của ~1^k + 2^k + ... + n^k = \sum_{i=1}^n i^k~ (sau khi chia lấy dư cho ~10^9 + 7~)

Link nộp bài: https://oj.vnoi.info/problem/bvd_tongmu

Hiển nhiên việc phương pháp duyệt với độ phức tạp ~O(nk)~ (hoặc ~O(n \log k)~ nếu dùng thuật toán tính lũy thừa nhanh) sẽ không giải quyết được bài toán này trong 1 giây, do số phép tính ~\approx nk \approx 10^{11} > 10^8~.

Ta thấy: Nếu ~k = 1~ thì ta có công thức tổng quát tính là ~\frac{n(n+1)}{2}~, nếu ~k = 2, 3~ thì trong sách giáo khoa Toán lớp 11 (cũ) cũng đã giới thiệu công thức tổng quát tính là một hàm bậc ~3~ và ~4~. Chính vì thế mà ta dự đoán được với mỗi ~k~ thì sẽ có một công thức tính dạng hàm số bậc ~k+1~ có biến ~n~.

Vậy ta mò công thức đó như thế nào?

Gọi công thức cần tìm là ~f(n)~

Đầu tiên, ta đã biết qua hai điểm không cùng hoành độ thì ta sẽ xác định được một đường thẳng (có phương trình là hàm bậc nhất), qua ba điểm không cùng hoành độ ta cũng xác định được một đường parabol (có phương trình là hàm bậc hai), vì vậy ta dự đoán được để tìm phương trình bậc ~(k+1)~ thì ta cần lấy ~(k+2)~ điểm không cùng hoành độ. Dễ nhất là ta lấy các điểm ~(1, f(1)), (2, f(2)), ..., (k+2, f(k+2))~. Mỗi điểm ta cần tối đa ~O(k \log k)~ phép tính để tính (dùng thuật toán tính lũy thừa nhanh), vì vậy ta sẽ mất ~O(k^2 \log k)~ phép tính cho bước này.

Ý tưởng chủ đạo của nội suy Lagrange để tìm một hàm đa thức đi qua các điểm là tách hàm ~f(n)~ thành tổng của các hàm ~g_1(n), g_2(n), ..., g_{k+2}(n)~ sao cho

~g_i(n) = \begin{cases} f(i) \text{ nếu } n = i \\ 0 \text{ nếu } n \text { là các hoành độ khác} \end{cases}~

Dễ thấy khi tìm được các hàm ~g~ thỏa mãn điều kiện trên thì ~f(n) = g_1(n) + g_2(n) + ... + g_{k+2}(n)~

Bây giờ ta tìm ~g_i(n)~ như thế nào.

Do ~g_i(n) = 0~ khi ~n~ là các hoành độ khác, điều này nghĩa là ~1, 2, 3, ..., i-1, i+1, ..., k+2~ (tất cả các số nguyên từ ~1~ đến ~k+2~ mà khác ~i~) là các nghiệm của ~g_i(n) = 0~. Do đó, ~g_i(n)~ sẽ có dạng ~\boxed{g_i(n) = (n-1)(n-2)(n-3)...(n-(i-1))(n-(i+1))...(n-(k+2))m_i}~ với ~m_i \in \mathbb{R}~ là một hằng số nào đó mà ta chưa biết. Ở đây có hai điểm đáng chú ý

• Cái này giống mẹo giải phương trình bằng cách nhẩm nghiệm ~x_0~, rồi cố gắng phân tích ra nhân tử ~(x - x_0)~ để ta có thể tính luôn nghiệm này và đơn giản hóa bài toán.
• Lí do ta biết ~m_i~ là hằng số là do phần ~(n-1)(n-2)(n-3)...~ là đã một đa thức bậc ~k+1~, mà ~f(n)~ là một đa thức bậc ~k+1~ nên ~m_i~ không được làm tăng bậc của ~g_i(n)~, hay nói cách khác, ~m_i~ là đa thức bậc 0 và là hằng số.

Để tính ~m_i~, ta sử dụng tính chất ~g_i(n) = f(i) \text{ nếu } n = i~.

Ta có: ~g_i(i) = f(i) \Rightarrow (i-1)(i-1)...(i-(i-1))(i - (i+1))...(i - (k+2))m_i = f(i) \Rightarrow \boxed{m_i = \frac{f(i)}{(i-1)(i-1)...(i-(i-1))(i - (i+1))...(i - (k+2))}}~

(ở phần mẫu không có thừa số nào bằng ~0~ vì hoành độ các điểm phân biệt, nên ta có thể chia được)

Vậy từ ~m_i~ ta ra được công thức của ~g_i(n)~, từ ~g_i(n)~ ta ra được công thức của ~f(n)~. Việc tính ~g_i(n)~ sẽ tốn ~O(k)~ phép tính, và ta cần tính ~k+2~ hàm như vậy, nên theo quy tắc nhân, độ phức tạp cuối cùng của bước này là ~O(k^2)~ (hoàn toàn không phụ thuộc vào ~n~)

Tổng độ phức tạp là ~O(k^2\log k + k^2) = O(\max(k^2\log k, k^2)) = O(k^2\log k)~

~k^2\log k \approx 10^7 < 10^8~, nên thuật toán sử dụng công thức ta mò ra bằng nội suy Lagrange có thể chạy được trong 1 giây.