矩陣乘法與斐波那契數列

前言

這篇文章屬於矩陣乘法的提高篇，雖然會對基礎知識進行講解，不過建議先進行學習后再來閱讀。
不保證能對您的水平帶來多大的提高，但一般來說會有的。

正文：

\(ps\):以下文章小寫字母及希臘字母代表一個實數，大寫字母代表矩陣，\(f_i\)代表斐波那契數列的第\(i\)項。

\(Part.1\) 矩陣運算

\(1\).加減法
若\(C=A\pm B\),則：

\[C_{i,j}=A_{i,j}\pm B_{i,j} \]

代碼實現：

struct jz
{
	long long a[105][105];
};
//兩個n*n的矩陣
jz add(jz x,jz y,int n)
{
	jz c;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) c.a[i][j]=0;	
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) c.a[i][j]=x.a[i][j]+y.a[i][j];
	}
	return c;
} 

\(2\).數乘
若\(C=kA\)，則：

\[C_{i,j}=kA_{i,j} \]

代碼實現：

struct jz
{
	long long a[105][105];
};
//一個n*n的矩陣
jz mathmul(jz x,long long k,int n)
{
	jz c;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) c.a[i][j]=0;	
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) c.a[i][j]=x.a[i][j]*k;
	}
	return c;
} 

\(3\).乘法
若\(C=A×B\)，則：

\[C_{i,j}=\sum\limits_{k=1}^n A_{i,k}×B_{k,j} \]

(這里的\(n\)代表\(A\)的行數或\(B\)的列數，當這兩個值不同時，相乘無意義)
代碼實現：

struct jz
{
	long long a[105][105];
};
//兩個n*n的矩陣
jz mul(jz x,jz y，int n)
{
	jz c;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) c.a[i][j]=0;	
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			for(int k=1;k<=n;k++)
			{
				c.a[i][j]=c.a[i][j]%MOD+x.a[i][k]*y.a[k][j]%MOD;
			}
		}
	}
	return c;
}

\(4\).求逆
不會，先咕着。

\(Part.2\) 矩陣快速冪

我們用自己定義的乘法跑快速冪即可。
例題題目鏈接：P3390 【模板】矩陣快速冪
代碼實現：

#include<iostream>
#include<cstdio>
using namespace std;
struct jz
{
	long a[105][105];
};
jz o;
const int MOD=1e9 +7;
int n;
jz yu;
long long  ci;
jz mul(jz x,jz y)
{
	jz c;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) c.a[i][j]=0;	
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			for(int k=1;k<=n;k++)
			{
				c.a[i][j]=c.a[i][j]%MOD+x.a[i][k]*y.a[k][j]%MOD;
			}
		}
	}
	return c;
}
jz quickpow(jz cc,long long k)
{
	jz ans=o,base=cc;
	while(k)
	{
		if(k&1) ans=mul(ans,base);
		base=mul(base,base);
		k>>=1;
	}
	return ans;
}
inline int read()
{
	int x=0,w=1;
	char c=getchar();
	while(c>'9'||c<'0')
	{
		if(c=='-')w=-1;
		c=getchar();
	}
	while(c<='9'&&c>='0')
	{
		x=(x<<3)+(x<<1)+c-'0';
		c=getchar();
	}
	return x*w;
}
int main()
{
   scanf("%d%lld",&n,&ci);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) yu.a[i][j]=read();
	}
	for(int i=1;i<=n;i++) o.a[i][i]=1;
	//這里的o是一個單位矩陣，大家可以想想為什么是這樣的
	jz d=quickpow(yu,ci);
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) 
		{
		    printf("%lld ",d.a[i][j]%MOD);
 		} 		
		printf("\n");
	}
	return 0;
} 

記得取模即可。

\(Part.3\) 疑惑

你可能會問：這是什么破算法，他能干啥？
剛學的我何嘗沒有這樣的疑問，不過現在了解之后，不會讓你們再有疑問了。
矩陣乘法是一種巧妙地方式將加法轉化成乘法的方式，以便在較短的方式解決遞推問題。
如果還是不太清楚，那我們就看些例子吧。

\(Part.4\) 簡單好啃的栗子

洛谷是有板子的：
題目鏈接:P1939 【模板】矩陣加速（數列）
對於這種加法形遞推式，一般都可以使用矩陣乘法加速遞推。
我們做矩陣乘法的原理是用一個矩陣封裝有用的變量，使用乘法實現多元加法，得到這個參數的下一項。
我們開始考慮要維護什么，顯然是\(a_i\)辣。
那我們考慮如何由\(a_{i-1}\)得到\(a_i\).
那我們先寫一個不全的遞推式：

\[a_{i-1}\longrightarrow a_i \]

然后嘗試在左邊加上某個（些）東西使等號成立。
很簡單：

\[a_{i-1}+a_{i-3}=a_i \]

顯然要考慮維護\(a_{i-3}\),並且記得去維護下一個下標的\(a_{i-3}\)，顯然是\(a_{i-2}\)。
接着考慮維護\(a_{i-2}\)的下一個是\(a_{i-1}\)，這個我們已經保存了，那么就好了。
我們嘗試用找到的規律畫一張表，標出相應的系數：

	\(a_{i-1}\)	\(a_{i-2}\)	\(a_{i-3}\)
\(a_{i-1}\longrightarrow a_i\)	\(1\)	\(0\)	\(1\)
\(a_{i-2}\longrightarrow a_{i-1}\)	\(1\)	\(0\)	\(0\)
\(a_{i-3}\longrightarrow a_{i-2}\)	\(0\)	\(1\)	\(0\)

這是我們構造和表格一樣的矩陣\(A\)：

\[\begin{pmatrix}1&0&1\\1&0&0\\0&1&0\end{pmatrix} \]

和：

\[\begin{pmatrix}a_{i-1}\\a_{i-2}\\a_{i-3}\end{pmatrix} \]

模擬下矩陣乘法，就會發現這就是我們所有想要的運算。
我們考慮從\(i=3\)開始，遞推\(k-2\)次，即把\(A\)進行\(k-2\)次冪。
然后乘上基礎的矩陣得到的值即可得到答案。
為防止文章過長，代碼自己寫吧，我就不給了。

\(Part.5\) 活學活用

再給你一道題：
題目鏈接：P1962 斐波那契數列
構造轉移矩陣：

\[\begin{pmatrix}1&1\\1&0\end{pmatrix} \]

太簡單了，再來一道：
題目鏈接;P1349 廣義斐波那契數列
轉移矩陣：

\[\begin{pmatrix}p&q\\1&0\end{pmatrix} \]

初始矩陣：

\[\begin{pmatrix}a_2\\a_1\end{pmatrix} \]

直接做就好了。

\(Part.6\) 嘗試搞事情

其實斐波那契數列有很多強勢的做法，這里只敘述兩種：
\(1\).我自己\(yy\)的。
比較久遠，大家可以去我的博客查看。-->快戳這兒，戳這兒
這里給出關鍵的式子，並命名為“斐波那契數列性質\(1\)”;

\[f_n=f_af_{n-a+1}+f_{a-1}f_{n-a} \]

可以用斐波那契數列定義（遞推式）去證明，不再贅述了。
\(2\).較為普及的算法。
我們稱之為“斐波那契數列性質\(2\)”：

\[\begin{cases}f_{2n}=f_{n+1}^2-f_{n-1}^2\\f_{2n+1}=f_n^2+f_{n+1}^2\end{cases} \]

第一個式子可以化簡：

\[\begin{aligned}f_{n+1}^2-f_{n-1}^2&=(f_{n+1}+f_{n-1})(f_{n+1}-f_{n-1})\cr &=(f_n+f_{n-1}+f_{n-1})f_n\cr&=(f_n+2f_{n-1})f_n\end{aligned} \]

數學歸納法證明。
首先驗證\(0<n<5\)時成立。
設第\(2n\)項之前的所有項都成立(不包含)。
那么:

\[\begin{aligned}f_{2n}&=f_{2n-2}+f_{2n-1}\cr &=f_n^2+f_{n-2}+f_n^2+f_{n-1}^2\cr&=f_n(f_{n-3}+2f_n)\cr&=f_n(f_{n-3}+f_n+f_{n-1}+f_{n-2})\cr&=f_n(f_n+2f_{n-1})\end{aligned} \]

得證。
另一種情況：

\[\begin{aligned}f_{2n+1}&=f_{2n}+f_{2n-1}\cr&=f_{n+1}^2-f_{n-1}^2+f_n^2+f_{n-1}^2\cr&=f_n^2+f_{n+1}^2\end{aligned} \]

得證。
其實很簡單的了，大頭在后面。

\(Part.7\) 斐波那契數列的通項公式

眾所周知：

\[f_n=\dfrac{1}{\sqrt{5}}((\dfrac{1+\sqrt{5}}{2})^n-(\dfrac{1-\sqrt{5}}{2})^n) \]

根據二項式定理展開可以得到（稱為斐波那契數列性質\(3\)）：

\[f_n=\dfrac{\sum\limits_{i=1}^n C^i_n×5^{\tfrac{n-1}{2}}×(i\%2)}{2^{n-1}} \]

暴力展開即可證明，不再贅述。
我們要用兩種方式證明通向公式：

方法一：生成函數法（不懂可以跳過）：

設斐波那契數列的生成函數為\(G(x)\)，即：

\[G(x)=\sum\limits_{i\geqslant 1}f_ix^i \]

然后套路的減一下：

\[G(x)-xG(x)=x+x^2G(x) \]

\(ps\)：上式靠手玩。

\[G(x)=\dfrac{x}{1-x-x^2} \]

因式分解：

\[G(x)=\dfrac{x}{(1-\dfrac{1-\sqrt{5}}{2}x)(1-\dfrac{1+\sqrt{5}}{2}x)} \]

博主只會暴算和亂湊，誰有好方法鴨？
然后裂項：

\[G(x)=-\dfrac{1}{\sqrt{5}}\dfrac{1}{(1-\dfrac{1-\sqrt{5}}{2}x)}+\dfrac{1}{\sqrt{5}}\dfrac{1}{(1-\dfrac{1+\sqrt{5}}{2}x)} \]

我再\(bb\)一下裂項吧：
對於上例,可以設\(a,b\),令：

\[a\dfrac{1}{(1-\dfrac{1-\sqrt{5}}{2}x)}+b\dfrac{1}{(1-\dfrac{1+\sqrt{5}}{2}x)}=G(x) \]

暴力解得\(a,b\)；
分治前面這一項，並把系數\(-\dfrac{1}{\sqrt{5}}\)提出來，要處理的就是這么一個玩意：

\[\dfrac{1}{(1-\dfrac{1-\sqrt{5}}{2}x)} \]

考慮到生成函數是等比數列求和的形式，把公比設成\(q\)，那么：

\[\dfrac{f_1(1-q^n)}{1-q}=\dfrac{1}{(1-\dfrac{1-\sqrt{5}}{2}x)} \]

把已知與極限代入可得：

\[\dfrac{1}{1-q}=\dfrac{1}{(1-\dfrac{1-\sqrt{5}}{2}x)} \]

解得：

\[q=\dfrac{1-\sqrt{5}}{2} \]

后面那些同理即得：

\[f_n=\dfrac{1}{\sqrt{5}}((\dfrac{1+\sqrt{5}}{2})^n-(\dfrac{1-\sqrt{5}}{2})^n) \]

這啥玩意，還這么難算？

方法二：特征根法：

首先明確所有線性遞推數列（學名“線性常系數齊次遞推”）都是等比數列。
我們有\(f_n=f_{n-1}+f_{n-2}\),那么有公比\(a\)滿足\(f_n=a^n\)，所有的\(a\)即為遞推式的特征根。
顯然：

\[a^n=a^{n-1}+a^{n-2} \]

同時除以\(a^{n-2}\)，得：

\[a^2-a-1=0 \]

\[a_1=\dfrac{1+\sqrt{5}}{2}\;,\;a_2=\dfrac{1-\sqrt{5}}{2} \]

那么有

\[f_n=xa_1^n+ya_2^n \]

把\(f_0=0,f_1=1\),分別帶入，得：

\[x=\dfrac{\sqrt{5}}{5}\;,\;y=-\dfrac{\sqrt{5}}{5} \]

同樣得到通項公式。

\(Part.8\) 斐波那契數列的幾何意義

給一張圖，奧妙無窮。

\(Part.9\) 稍有思維難度的矩陣題

上面扯出矩陣了，我們再扯回來。
注意：
\(1\).下面的題都是口胡題解，如果有鍋請積極提出，另外沒有代碼。
\(2\).我會“布置”一些作業，如果你覺得太屑了，就不做了吧。
\(3\).下面所有運算都在\(\bmod\)一個數的意義下，我就不寫了。
\(4\).數據范圍：\(n\leqslant 10^{18}\)。

\(T1\).

求\(\sum\limits_{i=1}^n f_i\) 。
做法一:
構造：

\[\begin{pmatrix}\sum\limits_{i=1}^{n-1}f_i\\f_{n-1}\\f_{n-2}\end{pmatrix}×\begin{pmatrix}1&1&1\\0&1&1\\0&1&0\end{pmatrix} \]

做法二:
利用斐波那契數列性質\(4\)：

\[\sum\limits_{i=1}^n f_i=f_{n+2}-1 \]

證明：

\[\sum\limits_{i=1}^n f_i=f_3+\sum\limits_{i=3}^n f_i=f_3+f_5+\sum\limits_{i=5}^n f_i=...... \]

分類討論：
當\(n\)為奇數時：

\[\sum\limits_{i=1}^n f_i=\sum\limits_{i=1}^{2i+1\leqslant n}f_{2i+1}+f_n \]

假設相等，那么：

\[\sum\limits_{i=1}^{2i+1\leqslant n}f_{2i+1}+f_n=f_{n+2}-1 \]

\[\sum\limits_{i=1}^{2i+1\leqslant n}f_{2i+1}+f_n=f_{n}+f_{n+1}-1 \]

\[\sum\limits_{i=1}^{2i+1\leqslant n}f_{2i+1}=f_{n+1}-1 \]

\[\sum\limits_{i=1}^{2i+1\leqslant n}f_{2i+1}=f_{n}+f_{n-1}-1 \]

我們會發現無限拆下去，出現了：

\[f_3=f_4-1 \]

的情況，並且這是成立的。
對於\(n\)為偶數的情況，也可以用同樣的方法驗證結論是成立的。
其實上面的證明都是吃多了的做法！
筆者在散步時想出了簡單易懂的證明方法。
考慮數學歸納法。
我們驗證\(n=1\)時成立。
那我們假設\(n\)及以前都成立：

\[\sum\limits_{i=1}^nf_i=f_{n+2}-1 \]

那么：

\[\sum\limits_{i=1}^{n+1}f_i=f_{n+2}+f_{n+1}-1=f_{n+3}-1 \]

得證。
然后用矩陣乘法求出\(f_{n+2}\)就好了。
作業\(1\):用幾何法證明利用斐波那契數列性質\(4\)。

\(T2\).

求\(\sum\limits_{i=1}^nf_i^2\)。
有難度，但是可以嘗試：

\[\sum\limits_{i=1}^{n-1}f_i^2\longrightarrow \sum\limits_{i=1}^{n}f_i^2 \]

顯然需要一個\(f_i^2\)。
考慮：

\[f_{n-1}^2\longrightarrow f_n^2 \]

由於：

\[\begin{aligned}f_n^2&=(f_{n-1}+f_{n-2})^2\cr&=f_{n-1}^2+f_{n-2}^2+2f_{n-1}f_{n-2}\end{aligned} \]

那么考慮維護\(f_{n-2}^2\)和\(f_{n-1}f_{n-2}\)。
由於\(f_{n-2}^2\)到\(f_{n-1}^2\)可以直接繼承，不用管了。
又因為：

\[f_{n}f_{n-1}=f_{n-1}(f_{n-1}+f_{n-2})=f_{n-1}^2+f_{n-1}f_{n-2} \]

都是可以繼承的，那么我們構造矩陣：

\[\begin{pmatrix}\sum\limits_{i=1}^{n-1}f_i^2\\f_{n-1}^2\\f_{n-2}^2\\f_{n-1}f_{n-2}\end{pmatrix}×\begin{pmatrix}1&1&1&2\\0&1&1&2\\0&1&0&0\\0&1&0&1\end{pmatrix} \]

好難啊！
我們可以證明斐波那契數列性質\(5\)來解決問題。
我們設\(S(n)=\sum\limits_{i=1}^nf_i^2,C(n)=\sum\limits_{i=3}^nf_{i-1}f_{i-2}\)
我們發現：

\[C(n)=\sum\limits_{i=3}^{n}f_{i-1}f_{i-2}=\sum\limits_{i=3}^{n}(f_{i-2}+f_{i-3})f_{i-2}=S(n-2)+C(n-1) \]

兩邊差分一下得：

\[S(n-2)=f_{n-1}f_{n-2} \]

很容易拓展到：

\[\sum\limits_{i=1}^nf_i^2=f_nf_{n+1} \]

算出\(f_n,f_{n+1}\)即可。
作業\(2\).用幾何法證明斐波那契數列性質\(5\)。

\(T3\).

求\(\sum\limits_{i=1}^n(n-i+1)f_i\)。
帶變化性系數，一眼很不可做。
但是發現他就是這么個東西：

\[\sum\limits_{i=1}^n \sum\limits_{j=1}^if_j \]

根據斐波那契數列性質\(4\):

\[\begin{aligned}\text{原式}&=\sum_{i=3}^{n+2} f_i-n\cr&=\sum\limits_{i=1}^{n+2}-n-2\cr&=f_{n+4}-n-3\end{aligned} \]

利用矩陣直接計算即可。

\(T4\).

把\(n\)個方塊排成一排，共有紅、藍、黃、綠四種顏色。求紅色與綠色方塊個數為偶數的方案數（不考慮順序，只考慮各種顏色的個數）。
一道好題。
方法一：
設當有\(i\)個方塊時，記紅綠都為偶數的方案數為\(a_i\),一個為偶數的方案數為\(b_i\),都不是的方案數為\(c_i\)。
容易得到：

\[\begin{cases}a_i=2a_{i-1}+b_{i-1}\\b_{i}=2a_{i-1}+2b_{i-1}+2c_{i-1}\\c_{i}=b_{i-1}+2c_{i-1}\end{cases} \]

那么構造：

\[\begin{pmatrix}a_{i-1}\\b_{i-1}\\c_{i-1}\end{pmatrix}×\begin{pmatrix}2&1&0\\2&2&2\\0&1&2\end{pmatrix} \]

方法二：
學過生成函數的同學都知道這是指數生成函數的入門題。
容易得到：

\[\begin{aligned}g^{(e)}&=(1+\dfrac{x^2}{2!}+\dfrac{x^4}{4!}+......)^2(1+\dfrac{x}{1!}+\dfrac{x^2}{2!}+......)^2\cr&=(\dfrac{1}{2}(e^x+e^{-x}))^2×e^{2x}\cr&=\dfrac{1}{4}(e^{4x}+2e^{2x}+1)\end{aligned} \]

容易得到通項公式：

\[h_i=\dfrac{4^i+2^{i+1}}{4}\;\;\;i>0 \]

emmm...直接快速冪就好了。

\(Part.10\) 另一類模型

題目鏈接：CF1182E Product Oriented Recurrence
帶上乘法就很難搞了，不過這是一種模型。
考慮維護系數,我們記(\(f\)不再是斐波那契數列)：

\[f_i=c^{w_i}f_1^{x_i}f_2^{y_i}f_3^{z_i} \]

容易得到遞推式：

\[\begin{cases}w_i=w_{i-1}+w_{i-2}+w_{i-3}+2i-4\\x_i=x_{i-1}+x_{i-2}+x_{i-3}\\y_{i}=y_{i-1}+y_{i-2}+y_{i-3}\\z_i=z_{i-1}+z_{i-2}+z_{i-3}\end{cases} \]

我們不妨讓數列從第\(4\)個開始編號，即：

\[f_i=c^{w_i-3}f_1^{x_i-3}f_2^{y_i-3}f_3^{z_i-3} \]

然后構造一大堆矩陣,發現\(x,y,z\)比較套路，就不說了。
對於\(w_i\)，構造：

\[\begin{pmatrix}w_{i-1}\\w_{i-2}\\w_{i-3}\\2i\\2\end{pmatrix}×\begin{pmatrix}1&1&1&1&0\\1&0&0&0&0\\0&1&0&0&0\\0&0&0&1&1\\0&0&0&0&1\end{pmatrix} \]

我們再運用拓展歐拉定理，來優化系數：不會的戳這里
由於\(\varphi(1e9+7)=1e9+6\),那我們把系數\(\bmod 1e9+6\)就好了。
但我們記得當數大於\(\varphi(m)\)時，要加\(\varphi(m)\)，我們打個表判斷就好，事實證明超過\(1e9+6\)的\(x_i,y_i,z_i,w_i\)對應的數列項數(從一開始編號，即把矩陣數列下標\(+1\).)分別為\(38,38,37,35\)，判斷一下即可。
要注意定義矩陣后清空，否則會死的很慘（隨機賦值我佛了）。
代碼實現：

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;

#define ll long long 
#define read(x) scanf("%lld",&x)
#define MOD 1000000006
#define MODD 1000000007

ll c,f1,f2,f3,rt;
ll ansx,ansy,ansz,answ;
ll ans=1;
struct mat
{	
	ll a[10][10];	
}e;

mat mul(mat x,mat y,int n,int mod)
{
	mat c;
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++) c.a[i][j]=0;	
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			for(int k=1;k<=n;k++)
			{
				c.a[i][j]=c.a[i][j]%mod+x.a[i][k]*y.a[k][j]%mod;
			}
		}
	}
	return c;
}

mat qp(mat cc,ll k,int n,int mod)
{
	mat ans=e,base=cc;
	while(k)
	{
		if(k&1) ans=mul(ans,base,n,mod);
		base=mul(base,base,n,mod);
		k>>=1;
	}
	return ans;
}

ll quickpow(ll a,ll b,ll mod)
{
	ll ans=1,base=a%mod;
	while(b)
	{
		if(b&1) ans=ans*base%mod;
		b>>=1;
		base=base*base%mod;	
	}	
	return ans%mod;
}

int main()
{
	read(rt),read(f1),read(f2),read(f3),read(c);
	//mod 1000000006;
	for(int i=1;i<=5;i++) for(int j=1;j<=5;j++) e.a[i][j]=0;
	for(int i=1;i<=5;i++) e.a[i][i]=1;
	mat x,xx;
	if(rt<=6) ansx=(rt<=5)?1:2;
	else
	{
		for(int i=1;i<=3;i++) for(int j=1;j<=3;j++) x.a[i][j]=0;
		for(int i=1;i<=3;i++) for(int j=1;j<=3;j++) xx.a[i][j]=0;
		x.a[1][1]=x.a[1][2]=x.a[1][3]=x.a[2][1]=x.a[3][2]=1;
		xx.a[1][1]=2,xx.a[2][1]=xx.a[3][1]=1;
		mat op1=qp(x,rt-6,3,MOD);
		op1=mul(op1,xx,3,MOD);
		ansx=op1.a[1][1]%MOD;
	}
	mat y,yy;
	if(rt<=6) ansy=rt-3;
	else 
	{
		for(int i=1;i<=3;i++) for(int j=1;j<=3;j++) y.a[i][j]=0;
		for(int i=1;i<=3;i++) for(int j=1;j<=3;j++) yy.a[i][j]=0;
		y.a[1][1]=y.a[1][2]=y.a[1][3]=y.a[2][1]=y.a[3][2]=1;
		yy.a[1][1]=3,yy.a[2][1]=2,yy.a[3][1]=1;
		mat op2=qp(y,rt-6,3,MOD);
		op2=mul(op2,yy,3,MOD);
		ansy=op2.a[1][1]%MOD;
	}
	mat z,zz;
	if(rt<=6) ansz=pow(2,rt-4);
	else
	{
		for(int i=1;i<=3;i++) for(int j=1;j<=3;j++) z.a[i][j]=0;
		for(int i=1;i<=3;i++) for(int j=1;j<=3;j++) zz.a[i][j]=0;
		z.a[1][1]=z.a[1][2]=z.a[1][3]=z.a[2][1]=z.a[3][2]=1;
		zz.a[1][1]=4,zz.a[2][1]=2;zz.a[3][1]=1;
		mat op3=qp(z,rt-6,3,MOD);
		op3=mul(op3,zz,3,MOD);
		ansz=op3.a[1][1]%MOD;
	}
	mat w,ww;
	if(rt<=6){if(rt==4) answ=2;else if(rt==5) answ=6;else answ=14;}
	else
	{
		for(int i=1;i<=5;i++) for(int j=1;j<=5;j++) w.a[i][j]=0;
		for(int i=1;i<=5;i++) for(int j=1;j<=5;j++) ww.a[i][j]=0;
		w.a[1][1]=w.a[1][2]=w.a[1][3]=w.a[1][4]=1;
		w.a[2][1]=w.a[3][2]=1;
		w.a[4][4]=w.a[4][5]=w.a[5][5]=1;
		ww.a[1][1]=14,ww.a[2][1]=6,ww.a[3][1]=2,ww.a[4][1]=8,ww.a[5][1]=2;
		mat op4=qp(w,rt-6,5,MOD);
		op4=mul(op4,ww,5,MOD);
		answ=op4.a[1][1]%MOD;
	}
	if(rt>=38) ans=ans*quickpow(f1,ansx+MOD,MODD)%MODD,ans=ans*quickpow(f2,ansy+MOD,MODD)%MODD;
	else ans=ans*quickpow(f1,ansx,MODD)%MODD,ans=ans*quickpow(f2,ansy,MODD)%MODD;
	if(rt>=37) ans=ans*quickpow(f3,ansz+MOD,MODD)%MODD;
	else ans=ans*quickpow(f3,ansz,MODD)%MODD;
	if(rt>=35) ans=ans*quickpow(c,answ+MOD,MODD)%MODD;
	else ans=ans*quickpow(c,answ,MODD)%MODD;
	printf("%lld\n",ans%MODD);
	return 0;
}

這種模型的原理：
利用已知的數列元素，把后面每一項表示出來，這是通常是這個元素冪的積的形式，通常使用拓展歐拉定理優化冪次。
核心思想是：同底數冪相乘，底數不變，指數相加。

拓展：
求\(\prod\limits_{i=1}^nf_i\)(這里的\(f_i\)指上面定義的數列）
我們考慮維護\(\sum x_i,\sum y_i,\sum z_i,\sum w_i\)
最后把\(x,y,z\)總和加上\(1\)即可。

\(Part.11\) 結語

首先感謝@Miraclys 大佬的審稿。
垃圾博主tlx同學耗費將近\(8\)個小時完成了這篇博客的構思\(+\)敘述。
他很累，當然他要寫的東西還有很多，好在他很快樂於與他人分享學習成果。
可悲的是\(AFO\)離他越來越近了。
以后像這么長的文章不知還有沒有，喜歡的話就給他點鼓勵吧（愣着干啥，點贊啊！）