矩阵乘法与斐波那契数列

前言

这篇文章属于矩阵乘法的提高篇，虽然会对基础知识进行讲解，不过建议先进行学习后再来阅读。
不保证能对您的水平带来多大的提高，但一般来说会有的。

正文：

\(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\)离他越来越近了。
以后像这么长的文章不知还有没有，喜欢的话就给他点鼓励吧（愣着干啥，点赞啊！）