alan23273850 的刷題筆記
d644: 壞脾氣小小皮
Tags: fast-exponentiation, fibonacci, matrix
出處：https://zerojudge.tw/ShowProblem?problemid=d644
提交：https://zerojudge.tw/Submissions?problemid=d644&account=allllllan123456
問題：定義數列 $a_1=1$、$a_2=1$，對於 $k\ge3$，$a_k = \begin{cases}a_{k-1}+a_{k-2} &\mbox{if }k \not\equiv 1\\a_{k-1}+a_{k-2}-1 & \mbox{if }k \equiv 1\end{cases}\pmod 3$，現在給你一個 int 的正整數 $n$，請求出 $a_n\pmod{100019}$。
解法：看到要找類費式數列的一般項數高達 int 上界的值，就要直接聯想到矩陣乘法搭配快速冪。那為了程式上實作的方便，讓 index 從 0 起算，定義 $b_k=a_{k+1}\ (\text{for}\ k\ge0)$，那麼矩陣乘法的一般式 $(t\ge0)$ 為：$\begin{bmatrix}b_{3(t+1)+0}\\b_{3(t+1)+1}\\b_{3(t+1)+2}\\-1\end{bmatrix}=\begin{bmatrix}0 & 1 & 1 & 1\\0 & 1 & 2 & 1\\0 & 2 & 3 & 2\\0 & 0 & 0 & 1\end{bmatrix}\begin{bmatrix}b_{3t+0}\\b_{3t+1}\\b_{3t+2}\\-1\end{bmatrix}$，數列三項三項一組再加上一個位移常量總共四項，矩陣的第一個列 $[0\ \ 1\ \ 1\ \ 1]$ 正處理了 $k\equiv0\pmod 3$ 的項，而第二列根據定義理論上要是它上方兩個列的總和，但是它上面只有一列，因為再更上面一列其實就對應到右側行向量的第三項，於是可以取而代之地直接對同一列的第三項加一，至於第三列才如剛剛所說的可以直接是第一列和第二列的相加，最後第四列必須維持位移常量才行。根據此遞迴式，易知一般項 $b_k$ 的求法便是 $\begin{bmatrix}0 & 1 & 1 & 1\\0 & 1 & 2 & 1\\0 & 2 & 3 & 2\\0 & 0 & 0 & 1\end{bmatrix}^{\displaystyle\lfloor k/3\rfloor}\begin{bmatrix}1\\1\\2\\-1\end{bmatrix}$ 之結果行向量的第 $k\pmod3$ 個元素。
1#include <bits/stdc++.h>
2using namespace std;
3
4struct matrix {
5 int e[4][4] = {{0,1,1,1},{0,1,2,1},{0,2,3,2},{0,0,0,1}};
6 matrix operator*(const matrix& m2) const {
7 matrix ans;
8 for (int i=0; i<4; i++)
9 for (int j=0; j<4; j++) {
10 long long temp = 0;
11 for (int k=0; k<4; k++)
12 temp += (long long)e[i][k] * m2.e[k][j];
13 ans.e[i][j] = temp % 100019;
14 }
15 return ans;
16 }
17};
18
19int main() { ios_base::sync_with_stdio(false); cin.tie(nullptr); // IO 優化
20 matrix mb; int ini2[4], ini[]={1,1,2,-1}, n; cin >> n; n--;
21 for (unsigned t=n/3; t; t>>=1) {
22 if (t & 1) {
23 memset(ini2, 0, sizeof ini2);
24 for (int i=0; i<4; i++) {
25 long long temp = 0;
26 for (int k=0; k<4; k++)
27 temp += (long long)mb.e[i][k] * ini[k];
28 ini2[i] = temp % 100019;
29 }
30 memcpy(ini, ini2, sizeof ini2);
31 }
32 mb = mb * mb;
33 }
34 cout << ini[n % 3] << '\n';
35 return 0;
36}