Nội suy
Giới thiệu
Nội suy là phương pháp suy ra giá trị dữ liệu mới trong một khoảng bằng cách dựa vào các điểm dữ liệu rời rạc đã biết. Nội suy thường dùng trong khớp hàm.
Ví dụ với các điểm dữ liệu:
\(x\) \(0\) \(1\) \(2\) \(3\) \(4\) \(5\) \(6\)
\(f(x)\) \(0\) \(0.8415\) \(0.9093\) \(0.1411\) \(-0.7568\) \(-0.9589\) \(-0.2794\)
Trong đó \(f(x)\) chưa biết, nội suy có thể ước lượng các điểm chưa biết bằng cách khớp \(f(x)\) theo một dạng nào đó.
Ví dụ, ta có thể dùng hàm tuyến tính từng đoạn:
Cách này gọi là nội suy tuyến tính .
Ta cũng có thể khớp bằng đa thức:
Cách này gọi là nội suy đa thức .
Dạng tổng quát của nội suy đa thức:
Nội suy đa thức
Với \(n+1\) điểm \((x_0,y_0),(x_1,y_1),\dots,(x_n,y_n)\) đã biết, tìm đa thức dạng \(f(x)=\sum_{i=0}^n a_ix^i\) thỏa
\[
f(x_i)=y_i,\qquad\forall i=0,1,\dots,n
\]
.
Sau đây giới thiệu hai cách trong nội suy đa thức: Lagrange và Newton. Có thể chứng minh hai cách cho cùng kết quả.
Nội suy Lagrange
Do cần dựng hàm \(f(x)\) đi qua các điểm \(P_1(x_1, y_1), P_2(x_2,y_2),\cdots,P_n(x_n,y_n)\) . Trước hết đặt hình chiếu của điểm thứ \(i\) lên trục \(x\) là \(P_i^{\prime}(x_i,0)\) .
Xét dựng \(n\) hàm \(f_1(x), f_2(x), \cdots, f_n(x)\) sao cho với hàm \(f_i(x)\) , đồ thị của nó đi qua \(\begin{cases}P_j^{\prime}(x_j,0),(j\neq i)\\P_i(x_i,y_i)\end{cases}\) , khi đó hàm cần tìm là \(f(x)=\sum\limits_{i=1}^nf_i(x)\) .
Ta đặt \(f_i(x)=a\cdot\prod_{j\neq i}(x-x_j)\) , thay \(P_i(x_i,y_i)\) vào suy ra \(a=\dfrac{y_i}{\prod_{j\neq i} (x_i-x_j)}\) , nên
\[
f_i(x)=y_i\cdot\dfrac{\prod_{j\neq i} (x-x_j)}{\prod_{j\neq i} (x_i-x_j)}=y_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}
\]
Suy ra dạng Lagrange:
\[
f(x)=\sum_{i=1}^ny_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j}
\]
Cài đặt ngây thơ có độ phức tạp \(O(n^2)\) , có thể tối ưu xuống \(O(n\log^2 n)\) , xem Nội suy đa thức nhanh .
Luogu P4781【模板】拉格朗日插值 Cho \(n\) điểm \((x_i,y_i)\) và \(k\) , với \(\forall i,j\) có \(i\neq j \iff x_i\neq x_j\) và \(f(x_i)\equiv y_i\pmod{998244353}\) và \(\deg(f(x)) < n\) (định nghĩa \(\deg(0)=-\infty\) ), tính \(f(k)\bmod{998244353}\) .
Lời giải
Ở bài này chỉ cần \(f(k)\) nên có thể thế trực tiếp \(k\) vào công thức; đôi khi cần nhiều lần tính hoặc thao tác phức tạp hơn thì phải tìm hệ số của \(f\) . Mã dưới đây cho cách tính hệ số.
\[
f(k)=\sum_{i=1}^{n}y_i\prod_{j\neq i }\frac{k-x_j}{x_i-x_j}
\]
Cần tính nghịch đảo. Nếu tính riêng tử và mẫu rồi nhân tử với nghịch đảo của mẫu, cộng dồn vào đáp án, nút thắt không còn ở nghịch đảo; độ phức tạp \(O(n^2)\) .
Vì làm việc trên modulo \(998244353\) , tạm coi thời gian tính nghịch đảo là hằng số.
Mã
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 #include <iostream>
#include <vector>
constexpr int MOD = 998244353 ;
using LL = long long ;
int inv ( int k ) {
int res = 1 ;
for ( int e = MOD - 2 ; e ; e /= 2 ) {
if ( e & 1 ) res = ( LL ) res * k % MOD ;
k = ( LL ) k * k % MOD ;
}
return res ;
}
// 返回 f 满足 f(x_i) = y_i
// 不考虑乘法逆元的时间，显然为 O(n^2)
std :: vector < int > lagrange_interpolation ( const std :: vector < int > & x ,
const std :: vector < int > & y ) {
const int n = x . size ();
std :: vector < int > M ( n + 1 ), px ( n , 1 ), f ( n );
M [ 0 ] = 1 ;
// 求出 M(x) = prod_(i=0..n-1)(x - x_i)
for ( int i = 0 ; i < n ; ++ i ) {
for ( int j = i ; j >= 0 ; -- j ) {
M [ j + 1 ] = ( M [ j ] + M [ j + 1 ]) % MOD ;
M [ j ] = ( LL ) M [ j ] * ( MOD - x [ i ]) % MOD ;
}
}
// 求出 px_i = prod_(j=0..n-1, j!=i) (x_i - x_j)
for ( int i = 0 ; i < n ; ++ i ) {
for ( int j = 0 ; j < n ; ++ j )
if ( i != j ) {
px [ i ] = ( LL ) px [ i ] * ( x [ i ] - x [ j ] + MOD ) % MOD ;
}
}
// 组合出 f(x) = sum_(i=0..n-1)(y_i / px_i)(M(x) / (x - x_i))
for ( int i = 0 ; i < n ; ++ i ) {
LL t = ( LL ) y [ i ] * inv ( px [ i ]) % MOD , k = M [ n ];
for ( int j = n - 1 ; j >= 0 ; -- j ) {
f [ j ] = ( f [ j ] + k * t ) % MOD ;
k = ( M [ j ] + k * x [ i ]) % MOD ;
}
}
return f ;
}
int main () {
std :: ios :: sync_with_stdio ( false );
std :: cin . tie ( nullptr );
int n , k ;
std :: cin >> n >> k ;
std :: vector < int > x ( n ), y ( n );
for ( int i = 0 ; i < n ; ++ i ) std :: cin >> x [ i ] >> y [ i ];
const auto f = lagrange_interpolation ( x , y );
int v = 0 ;
for ( int i = n - 1 ; i >= 0 ; -- i ) v = (( LL ) v * k + f [ i ]) % MOD ;
std :: cout << v << '\n' ;
return 0 ;
}
Nội suy Lagrange với hoành độ là các số nguyên liên tiếp
Nếu biết hoành độ là các số nguyên liên tiếp, có thể nội suy \(O(n)\) .
Gọi đa thức cần tìm là \(f(x)\) , biết \(f(1),\cdots,f(n+1)\) (\(1\le i\le n+1\) ), xét công thức:
\[
\begin{aligned}
f(x)&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-x_j}{x_i-x_j}\\
&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-j}{i-j}
\end{aligned}
\]
Tích phía sau xét tử và mẫu riêng; tử là:
\[
\dfrac{\prod\limits_{j=1}^{n+1}(x-j)}{x-i}
\]
Mẫu tích \(i-j\) tách thành hai giai đoạn theo giai thừa:
\[
(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!
\]
Suy ra công thức nội suy khi hoành độ là \(1,\cdots,n+1\) :
\[
f(x)=\sum\limits_{i=1}^{n+1}(-1)^{n+1-i}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(i-1)!(n+1-i)!(x-i)}
\]
Tiền xử lý tích tiền/sau của \((x-i)\) , giai thừa và nghịch đảo giai thừa, rồi thay vào công thức, độ phức tạp \(O(n)\) .
Ví dụ CF622F The Sum of the k-th Powers
Cho \(n,k\) , tính \(\sum\limits_{i=1}^ni^k\) modulo \(10^9+7\) .
Lời giải
Đáp án là đa thức bậc \(k+1\) , nên có thể sàng tuyến tính \(1^i,\cdots,(k+2)^i\) rồi nội suy \(O(n)\) .
Cũng có thể từ kiến thức tổ hợp và sai phân suy ra:
\[
f(x)=\sum_{i=1}^{n+1}\binom{x-1}{i-1}\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}y_{j}=\sum\limits_{i=1}^{n+1}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(x-i)\cdot(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!}
\]
Mã
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 // By: Luogu@rui_er(122461)
#include <iostream>
using namespace std ;
constexpr int N = 1e6 + 5 , mod = 1e9 + 7 ;
int n , k , tab [ N ], p [ N ], pcnt , f [ N ], pre [ N ], suf [ N ], fac [ N ], inv [ N ], ans ;
int qpow ( int x , int y ) {
int ans = 1 ;
for (; y ; y >>= 1 , x = 1L L * x * x % mod )
if ( y & 1 ) ans = 1L L * ans * x % mod ;
return ans ;
}
void sieve ( int lim ) {
f [ 1 ] = 1 ;
for ( int i = 2 ; i <= lim ; i ++ ) {
if ( ! tab [ i ]) {
p [ ++ pcnt ] = i ;
f [ i ] = qpow ( i , k );
}
for ( int j = 1 ; j <= pcnt && 1L L * i * p [ j ] <= lim ; j ++ ) {
tab [ i * p [ j ]] = 1 ;
f [ i * p [ j ]] = 1L L * f [ i ] * f [ p [ j ]] % mod ;
if ( ! ( i % p [ j ])) break ;
}
}
for ( int i = 2 ; i <= lim ; i ++ ) f [ i ] = ( f [ i - 1 ] + f [ i ]) % mod ;
}
int main () {
cin . tie ( nullptr ) -> sync_with_stdio ( false );
cin >> n >> k ;
sieve ( k + 2 );
if ( n <= k + 2 ) return cout << f [ n ], 0 ;
pre [ 0 ] = suf [ k + 3 ] = 1 ;
for ( int i = 1 ; i <= k + 2 ; i ++ ) pre [ i ] = 1L L * pre [ i - 1 ] * ( n - i ) % mod ;
for ( int i = k + 2 ; i >= 1 ; i -- ) suf [ i ] = 1L L * suf [ i + 1 ] * ( n - i ) % mod ;
fac [ 0 ] = inv [ 0 ] = fac [ 1 ] = inv [ 1 ] = 1 ;
for ( int i = 2 ; i <= k + 2 ; i ++ ) {
fac [ i ] = 1L L * fac [ i - 1 ] * i % mod ;
inv [ i ] = 1L L * ( mod - mod / i ) * inv [ mod % i ] % mod ;
}
for ( int i = 2 ; i <= k + 2 ; i ++ ) inv [ i ] = 1L L * inv [ i - 1 ] * inv [ i ] % mod ;
for ( int i = 1 ; i <= k + 2 ; i ++ ) {
int P = 1L L * pre [ i - 1 ] * suf [ i + 1 ] % mod ;
int Q = 1L L * inv [ i - 1 ] * inv [ k + 2 - i ] % mod ;
int mul = (( k + 2 - i ) & 1 ) ? -1 : 1 ;
ans = ( ans + 1L L * ( Q * mul + mod ) % mod * P % mod * f [ i ] % mod ) % mod ;
}
cout << ans << '\n' ;
return 0 ;
}
Nội suy Newton
Nội suy Newton dựa trên sai phân bậc cao, ưu điểm là hỗ trợ thêm điểm mới trong \(O(n)\) .
Để chèn điểm trong \(O(n)\) , đặt:
\[
f(x)=\sum_{j=0}^n a_jn_j(x)
\]
với \(n_j(x):=\prod_{i=0}^{j-1}(x-x_i)\) gọi là cơ sở Newton (Newton basis).
Giải \(a_j\) sẽ được đa thức nội suy. Ta định nghĩa thương sai tiến (forward divided differences):
\[
\begin{aligned}
\lbrack y_k\rbrack & := y_k, & k=0,\dots,n, \\
[y_k,\dots,y_{k+j}] & := \dfrac{[y_{k+1},\dots,y_{k+j}]-[y_k,\dots,y_{k+j-1}]}{x_{k+j}-x_k}, & k=0,\dots,n-j,~j=1,\dots,n.
\end{aligned}
\]
Khi đó:
\[
\begin{aligned}
f(x)&=[y_0]+[y_0,y_1](x-x_0)+\dots+[y_0,\dots,y_n](x-x_0)\dots(x-x_{n-1})\\
&=\sum_{j=0}^n [y_0,\dots,y_j]n_j(x)
\end{aligned}
\]
Đây là dạng Newton. Cài đặt ngây thơ \(O(n^2)\) .
Nếu các điểm mẫu cách đều (tức \(x_i=x_0+ih\) , \(i=1,\dots,n\) ), ta có:
\[
[y_k,\dots,y_{k+j}]=\frac{1}{j!h^j}\Delta^{(j)}y_k,
\]
trong đó \(\Delta^{(j)}y_k\) là sai phân tiến (forward differences), định nghĩa:
\[
\begin{aligned}
\Delta^{(0)}y_k & := y_k, & k=0,\dots,n, \\
\Delta^{(j)}y_k & := \Delta^{(j-1)} y_{k+1}-\Delta^{(j-1)} y_k, & k=0,\dots,n-j,~j=1,\dots,n.
\end{aligned}
\]
Đặt \(x=x_0+sh\) , công thức Newton thành:
\[
f(x)=\sum_{j=0}^n \binom{s}{j}j!h^j[y_0,\dots,y_j]=\sum_{j=0}^n \binom{s}{j}\Delta^{(j)}y_0.
\]
Code (Luogu P4781【模板】拉格朗日插值 )
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
106
107
108
109
110
111 #include <cstdint>
#include <iostream>
#include <vector>
using namespace std ;
constexpr uint32_t MOD = 998244353 ;
struct mint {
uint32_t v_ ;
mint () : v_ ( 0 ) {}
mint ( int64_t v ) {
int64_t x = v % ( int64_t ) MOD ;
v_ = ( uint32_t )( x + ( x < 0 ? MOD : 0 ));
}
friend mint inv ( mint const & x ) {
int64_t a = x . v_ , b = MOD ;
if (( a %= b ) == 0 ) return 0 ;
int64_t s = b , m0 = 0 ;
for ( int64_t q = 0 , _ = 0 , m1 = 1 ; a ;) {
_ = s - a * ( q = s / a );
s = a ;
a = _ ;
_ = m0 - m1 * q ;
m0 = m1 ;
m1 = _ ;
}
return m0 ;
}
mint & operator += ( mint const & r ) {
if (( v_ += r . v_ ) >= MOD ) v_ -= MOD ;
return * this ;
}
mint & operator -= ( mint const & r ) {
if (( v_ -= r . v_ ) >= MOD ) v_ += MOD ;
return * this ;
}
mint & operator *= ( mint const & r ) {
v_ = ( uint32_t )(( uint64_t ) v_ * r . v_ % MOD );
return * this ;
}
mint & operator /= ( mint const & r ) { return * this = * this * inv ( r ); }
friend mint operator + ( mint l , mint const & r ) { return l += r ; }
friend mint operator - ( mint l , mint const & r ) { return l -= r ; }
friend mint operator * ( mint l , mint const & r ) { return l *= r ; }
friend mint operator / ( mint l , mint const & r ) { return l /= r ; }
};
template < class T >
struct NewtonInterp {
// {(x_0,y_0),...,(x_{n-1},y_{n-1})}
vector < pair < T , T >> p ;
// dy[r][l] = [y_l,...,y_r]
vector < vector < T >> dy ;
// (x-x_0)...(x-x_{n-1})
vector < T > base ;
// [y_0]+...+[y_0,y_1,...,y_n](x-x_0)...(x-x_{n-1})
vector < T > poly ;
void insert ( T const & x , T const & y ) {
p . emplace_back ( x , y );
size_t n = p . size ();
if ( n == 1 ) {
base . push_back ( 1 );
} else {
size_t m = base . size ();
base . push_back ( 0 );
for ( size_t i = m ; i ; -- i ) base [ i ] = base [ i - 1 ];
base [ 0 ] = 0 ;
for ( size_t i = 0 ; i < m ; ++ i )
base [ i ] = base [ i ] - p [ n - 2 ]. first * base [ i + 1 ];
}
dy . emplace_back ( p . size ());
dy [ n - 1 ][ n - 1 ] = y ;
if ( n > 1 ) {
for ( size_t i = n - 2 ; ~ i ; -- i )
dy [ n - 1 ][ i ] =
( dy [ n - 2 ][ i ] - dy [ n - 1 ][ i + 1 ]) / ( p [ i ]. first - p [ n - 1 ]. first );
}
poly . push_back ( 0 );
for ( size_t i = 0 ; i < n ; ++ i ) poly [ i ] = poly [ i ] + dy [ n - 1 ][ 0 ] * base [ i ];
}
T eval ( T const & x ) {
T ans {};
for ( auto it = poly . rbegin (); it != poly . rend (); ++ it ) ans = ans * x + * it ;
return ans ;
}
};
int main () {
NewtonInterp < mint > ip ;
int n , k ;
cin >> n >> k ;
for ( int i = 1 , x , y ; i <= n ; ++ i ) {
cin >> x >> y ;
ip . insert ( x , y );
}
cout << ip . eval ( k ). v_ ;
return 0 ;
}
Nội suy Newton với hoành độ là các số nguyên liên tiếp
Ví dụ: tìm hệ số của đa thức \(f(x)=\sum_{i=0}^{3} a_ix^i\) , biết \(f(1)\) đến \(f(6)\) lần lượt là \(1, 5, 14, 30, 55, 91\) .
\[
\begin{array}{cccccccccccc}
1 & & 5 & & 14 & & 30 & & 55 & & 91 & \\
& 4 & & 9 & & 16 & & 25 & & 36 & \\
& & 5 & & 7 & & 9 & & 11 & \\
& & & 2 & & 2 & & 2 & \\
\end{array}
\]
Hàng đầu là các giá trị liên tiếp của \(f(x)\) ; mỗi hàng sau là hiệu giữa hai phần tử kề nhau của hàng trước. Quan sát: nếu thao tác đủ nhiều lần (với \(f(x)\) là đa thức), cuối cùng sẽ thu được một giá trị hằng.
Phần tử đầu của sai phân bậc \(i-1\) là \(\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j)\) , và đóng góp của nó vào \(f(k)\) là \(\binom{k-1}{i-1}\) lần.
\[
f(k)=\sum_{i=1}^n\binom{k-1}{i-1}\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j)
\]
Độ phức tạp \(O(n^2)\) .
Cài đặt trong C++
Từ C++20, thư viện chuẩn có std::midpointstd::lerp
Bài tập
Tài liệu tham khảo
Interpolation - Wikipedia Newton polynomial - Wikipedia
Update history Edit this page on GitHub! AtomAlpaca, billchenchina, caibyte, Chrogeek, Early0v0, EndlessCheng, Enter-tainer, Henry-ZHR, hly1204, hsfzLZH1, Ir1d, Ghastlcon, kenlig, Marcythm, megakite, Peanut-Tang, qwqAutomaton, qz-cqy, StudyingFather, swift-zym, swiftqwq, Tiphereth-A, TrisolarisHD, Watersail2005, x4Cx58x54, Xeonacid, xiaopangfeiyu, YanWQ-monad CC BY-SA 4.0 and SATA
