矩阵乘法与斐波那契数列

前言

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

正文：

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