This is the first activity of Computerization algorithm team. We introduced the method to find the n n n th term of the Fibonacci sequence, which mainly uses matrix exponentiation.
Problemโ
The Fibonacci sequence:
F n = { 0 , n = 0 1 , n = 1 F n โ 2 + F n โ 1 , n > 1 F_{n}=
\begin{cases}
0,&n=0\\
1,&n=1\\
F_{n-2}+F_{n-1},&n>1
\end{cases} F n โ = โฉ โจ โง โ 0 , 1 , F n โ 2 โ + F n โ 1 โ , โ n = 0 n = 1 n > 1 โ Given n n n , find F n ย modย 1 0 9 + 7 F_{n}\text{ mod }10^9+7 F n โ ย modย 1 0 9 + 7 ใ
Input constraints Memory limit Execution time 0 โค n โค 1 0 19 0\le n\le 10^{19} 0 โค n โค 1 0 19 64MB 1.0s
Solutionโ
The input size of 1 0 19 10^{19} 1 0 19 obviously prohibits any attempt to solve it with loops. Is there a better way than a simple O ( n ) \mathcal{O}(n) O ( n ) ? It turns out that with matrix exponentiation , we can achieve O ( log โก n ) \mathcal{O}(\log n) O ( log n ) . We observe that:
( F n + 1 F n + 2 ) = ( F n + 1 F n + F n + 1 ) = ( 0 1 1 1 ) ( F n F n + 1 ) \begin{pmatrix}F_{n+1}\\F_{n+2}\end{pmatrix}=\begin{pmatrix}F_{n+1}\\F_{n}+F_{n+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix} ( F n + 1 โ F n + 2 โ โ ) = ( F n + 1 โ F n โ + F n + 1 โ โ ) = ( 0 1 โ 1 1 โ ) ( F n โ F n + 1 โ โ )
This step is applicable to all recursive sequences, so it should be easily reached for an experienced candidate. Generally, for F n + 2 = a F n + b F n + 1 F_{n+2}=aF_{n}+bF_{n+1} F n + 2 โ = a F n โ + b F n + 1 โ , we have
( F n + 1 F n + 2 ) = ( F n + 1 a F n + b F n + 1 ) = ( 0 1 a b ) ( F n F n + 1 ) \begin{pmatrix}F_{n+1}\\F_{n+2}\end{pmatrix}=\begin{pmatrix}F_{n+1}\\aF_{n}+bF_{n+1}\end{pmatrix}=\begin{pmatrix}0&1\\a&b\end{pmatrix}\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix} ( F n + 1 โ F n + 2 โ โ ) = ( F n + 1 โ a F n โ + b F n + 1 โ โ ) = ( 0 a โ 1 b โ ) ( F n โ F n + 1 โ โ )
From the recursive definition,
( F n + m F n + m + 1 ) = ( 0 1 1 1 ) m ( F n F n + 1 ) \begin{pmatrix}F_{n+m}\\F_{n+m+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^m\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix} ( F n + m โ F n + m + 1 โ โ ) = ( 0 1 โ 1 1 โ ) m ( F n โ F n + 1 โ โ )
Substituting n = 0 n=0 n = 0 ,
( F m F m + 1 ) = ( 0 1 1 1 ) m ( F 0 F 1 ) \begin{pmatrix}F_{m}\\F_{m+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^m\begin{pmatrix}F_0\\F_1\end{pmatrix} ( F m โ F m + 1 โ โ ) = ( 0 1 โ 1 1 โ ) m ( F 0 โ F 1 โ โ )
Now the problem is transformed into finding the matrix raised to the m m m th power. If m = 2 0 a 0 + 2 1 a 1 + 2 2 a 2 + โฆ m=2^0a_0+2^1a_1+2^2a_2+\dots m = 2 0 a 0 โ + 2 1 a 1 โ + 2 2 a 2 โ + โฆ (representation in binary), then
( 0 1 1 1 ) m = ( ( 0 1 1 1 ) 1 ) a 0 ร ( ( 0 1 1 1 ) 2 ) a 1 ร ( ( 0 1 1 1 ) 4 ) a 2 โฆ \begin{pmatrix}0&1\\1&1\end{pmatrix}^m=\left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{1}\right)^{a_0}\times \left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{2}\right)^{a_1}\times \left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{4}\right)^{a_2}\dots ( 0 1 โ 1 1 โ ) m = ( ( 0 1 โ 1 1 โ ) 1 ) a 0 โ ร ( ( 0 1 โ 1 1 โ ) 2 ) a 1 โ ร ( ( 0 1 โ 1 1 โ ) 4 ) a 2 โ โฆ
The 2 k 2^k 2 k th powers of the original matrix can, in fact, be preprocessed. When m < 1 0 19 m<10^{19} m < 1 0 19 , k < log โก 2 1 0 19 < 64 k<\log_2 10^{19}<64 k < log 2 โ 1 0 19 < 64 , so we only need to store at most 63 matrices. In addition,
( 0 1 1 1 ) 2 k = ( 0 1 1 1 ) 2 k โ 1 ร ( 0 1 1 1 ) 2 k โ 1 \begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^k}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^{k-1}}\times \begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^{k-1}} ( 0 1 โ 1 1 โ ) 2 k = ( 0 1 โ 1 1 โ ) 2 k โ 1 ร ( 0 1 โ 1 1 โ ) 2 k โ 1
which implies that the powers can be attained within O ( log โก m ) \mathcal{O}(\log m) O ( log m ) time. This is the idea of fast matrix exponentiation: compute all 2 k 2^k 2 k th powers, and put those needed together.
Programโ
Below is the C++ code, where the most intractable part is probably implementation of matrix multiplication:
# include <iostream> # include <cmath> using namespace std ; struct mat { unsigned long long a [ 4 ] ; mat operator * ( mat o ) { mat t ; t . a [ 0 ] = ( this -> a [ 0 ] * o . a [ 0 ] + this -> a [ 1 ] * o . a [ 2 ] ) % 1000000007 ; t . a [ 1 ] = ( this -> a [ 0 ] * o . a [ 1 ] + this -> a [ 1 ] * o . a [ 3 ] ) % 1000000007 ; t . a [ 2 ] = ( this -> a [ 2 ] * o . a [ 0 ] + this -> a [ 3 ] * o . a [ 2 ] ) % 1000000007 ; t . a [ 3 ] = ( this -> a [ 2 ] * o . a [ 1 ] + this -> a [ 3 ] * o . a [ 3 ] ) % 1000000007 ; return t ; } } ; mat mat_pow [ 64 ] ; int fib ( unsigned long long k ) { mat tmp ; tmp . a [ 0 ] = 1 ; tmp . a [ 1 ] = 0 ; tmp . a [ 2 ] = 0 ; tmp . a [ 3 ] = 1 ; for ( int i = 0 ; i < 64 ; i ++ ) { if ( k & ( 1ull << i ) ) { tmp = tmp * mat_pow [ i ] ; } } return tmp . a [ 1 ] ; } int main ( ) { mat_pow [ 0 ] . a [ 0 ] = 0 ; mat_pow [ 0 ] . a [ 1 ] = 1 ; mat_pow [ 0 ] . a [ 2 ] = 1 ; mat_pow [ 0 ] . a [ 3 ] = 1 ; for ( int i = 1 ; i < 64 ; i ++ ) { mat_pow [ i ] = mat_pow [ i - 1 ] * mat_pow [ i - 1 ] ; } unsigned long long n ; cin >> n ; cout << fib ( n ) << endl ; return 0 ; }
The formulae for matrix multiplication are:
( a 0 a 1 a 2 a 3 ) ร ( b 0 b 1 b 2 b 3 ) = ( a 0 b 0 + a 1 b 2 a 0 b 1 + a 1 b 3 a 2 b 0 + a 3 b 2 a 2 b 1 + a 3 b 3 ) \begin{pmatrix}a_0&a_1\\a_2&a_3\end{pmatrix}\times\begin{pmatrix}b_0&b_1\\b_2&b_3\end{pmatrix}=\begin{pmatrix}a_0b_0+a_1b_2&a_0b_1+a_1b_3\\a_2b_0+a_3b_2&a_2b_1+a_3b_3\end{pmatrix} ( a 0 โ a 2 โ โ a 1 โ a 3 โ โ ) ร ( b 0 โ b 2 โ โ b 1 โ b 3 โ โ ) = ( a 0 โ b 0 โ + a 1 โ b 2 โ a 2 โ b 0 โ + a 3 โ b 2 โ โ a 0 โ b 1 โ + a 1 โ b 3 โ a 2 โ b 1 โ + a 3 โ b 3 โ โ )
( a 0 a 1 a 2 a 3 ) ร ( b 0 b 1 ) = ( a 0 b 0 + a 1 b 1 a 2 b 0 + a 3 b 1 ) \begin{pmatrix}a_0&a_1\\a_2&a_3\end{pmatrix}\times\begin{pmatrix}b_0\\b_1\end{pmatrix}=\begin{pmatrix}a_0b_0+a_1b_1\\a_2b_0+a_3b_1\end{pmatrix} ( a 0 โ a 2 โ โ a 1 โ a 3 โ โ ) ร ( b 0 โ b 1 โ โ ) = ( a 0 โ b 0 โ + a 1 โ b 1 โ a 2 โ b 0 โ + a 3 โ b 1 โ โ )