矩阵乘法与斐波那契数列

前言

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

正文：

(ps):以下文章小写字母及希腊字母代表一个实数，大写字母代表矩阵，(f_i)代表斐波那契数列的第(i)项。

(Part.1) 矩阵运算

(1).加减法
若(C=Apm 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}=sumlimits_{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("
");
	}
	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)：

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

和：

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

模拟下矩阵乘法，就会发现这就是我们所有想要的运算。
我们考虑从(i=3)开始，递推(k-2)次，即把(A)进行(k-2)次幂。
然后乘上基础的矩阵得到的值即可得到答案。
为防止文章过长，代码自己写吧，我就不给了。

(Part.5) 活学活用

再给你一道题：
题目链接：P1962 斐波那契数列
构造转移矩阵：

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

太简单了，再来一道：
题目链接;P1349 广义斐波那契数列
转移矩阵：

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

初始矩阵：

[egin{pmatrix}a_2\a_1end{pmatrix} ]

直接做就好了。

(Part.6) 尝试搞事情

其实斐波那契数列有很多强势的做法，这里只叙述两种：
(1).我自己(yy)的。
比较久远，大家可以去我的博客查看。-->快戳这儿，戳这儿
这里给出关键的式子，并命名为“斐波那契数列性质(1)”;

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

可以用斐波那契数列定义（递推式）去证明，不再赘述了。
(2).较为普及的算法。
我们称之为“斐波那契数列性质(2)”：

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

第一个式子可以化简：

[egin{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_ncr&=(f_n+2f_{n-1})f_nend{aligned} ]

数学归纳法证明。
首先验证(0<n<5)时成立。
设第(2n)项之前的所有项都成立(不包含)。
那么:

[egin{aligned}f_{2n}&=f_{2n-2}+f_{2n-1}cr &=f_n^2+f_{n-2}+f_n^2+f_{n-1}^2cr&=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} ]

得证。
另一种情况：

[egin{aligned}f_{2n+1}&=f_{2n}+f_{2n-1}cr&=f_{n+1}^2-f_{n-1}^2+f_n^2+f_{n-1}^2cr&=f_n^2+f_{n+1}^2end{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{sumlimits_{i=1}^n C^i_n×5^{ frac{n-1}{2}}×(i\%2)}{2^{n-1}} ]

暴力展开即可证明，不再赘述。
我们要用两种方式证明通向公式：

方法一：生成函数法（不懂可以跳过）：

设斐波那契数列的生成函数为(G(x))，即：

[G(x)=sumlimits_{igeqslant 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),令：

[adfrac{1}{(1-dfrac{1-sqrt{5}}{2}x)}+bdfrac{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).下面所有运算都在(mod)一个数的意义下，我就不写了。
(4).数据范围：(nleqslant 10^{18})。

(T1).

求(sumlimits_{i=1}^n f_i) 。
做法一:
构造：

[egin{pmatrix}sumlimits_{i=1}^{n-1}f_i\f_{n-1}\f_{n-2}end{pmatrix}×egin{pmatrix}1&1&1\0&1&1\0&1&0end{pmatrix} ]

做法二:
利用斐波那契数列性质(4)：

[sumlimits_{i=1}^n f_i=f_{n+2}-1 ]

证明：

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

分类讨论：
当(n)为奇数时：

[sumlimits_{i=1}^n f_i=sumlimits_{i=1}^{2i+1leqslant n}f_{2i+1}+f_n ]

假设相等，那么：

[sumlimits_{i=1}^{2i+1leqslant n}f_{2i+1}+f_n=f_{n+2}-1 ]

[sumlimits_{i=1}^{2i+1leqslant n}f_{2i+1}+f_n=f_{n}+f_{n+1}-1 ]

[sumlimits_{i=1}^{2i+1leqslant n}f_{2i+1}=f_{n+1}-1 ]

[sumlimits_{i=1}^{2i+1leqslant n}f_{2i+1}=f_{n}+f_{n-1}-1 ]

我们会发现无限拆下去，出现了：

[f_3=f_4-1 ]

的情况，并且这是成立的。
对于(n)为偶数的情况，也可以用同样的方法验证结论是成立的。
其实上面的证明都是吃多了的做法！
笔者在散步时想出了简单易懂的证明方法。
考虑数学归纳法。
我们验证(n=1)时成立。
那我们假设(n)及以前都成立：

[sumlimits_{i=1}^nf_i=f_{n+2}-1 ]

那么：

[sumlimits_{i=1}^{n+1}f_i=f_{n+2}+f_{n+1}-1=f_{n+3}-1 ]

得证。
然后用矩阵乘法求出(f_{n+2})就好了。
作业(1):用几何法证明利用斐波那契数列性质(4)。

(T2).

求(sumlimits_{i=1}^nf_i^2)。
有难度，但是可以尝试：

[sumlimits_{i=1}^{n-1}f_i^2longrightarrow sumlimits_{i=1}^{n}f_i^2 ]

显然需要一个(f_i^2)。
考虑：

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

由于：

[egin{aligned}f_n^2&=(f_{n-1}+f_{n-2})^2cr&=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} ]

都是可以继承的，那么我们构造矩阵：

[egin{pmatrix}sumlimits_{i=1}^{n-1}f_i^2\f_{n-1}^2\f_{n-2}^2\f_{n-1}f_{n-2}end{pmatrix}×egin{pmatrix}1&1&1&2\0&1&1&2\0&1&0&0\0&1&0&1end{pmatrix} ]

好难啊！
我们可以证明斐波那契数列性质(5)来解决问题。
我们设(S(n)=sumlimits_{i=1}^nf_i^2,C(n)=sumlimits_{i=3}^nf_{i-1}f_{i-2})
我们发现：

[C(n)=sumlimits_{i=3}^{n}f_{i-1}f_{i-2}=sumlimits_{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} ]

很容易拓展到：

[sumlimits_{i=1}^nf_i^2=f_nf_{n+1} ]

算出(f_n,f_{n+1})即可。
作业(2).用几何法证明斐波那契数列性质(5)。

(T3).

求(sumlimits_{i=1}^n(n-i+1)f_i)。
带变化性系数，一眼很不可做。
但是发现他就是这么个东西：

[sumlimits_{i=1}^n sumlimits_{j=1}^if_j ]

根据斐波那契数列性质(4):

[egin{aligned} ext{原式}&=sum_{i=3}^{n+2} f_i-ncr&=sumlimits_{i=1}^{n+2}-n-2cr&=f_{n+4}-n-3end{aligned} ]

利用矩阵直接计算即可。

(T4).

把(n)个方块排成一排，共有红、蓝、黄、绿四种颜色。求红色与绿色方块个数为偶数的方案数（不考虑顺序，只考虑各种颜色的个数）。
一道好题。
方法一：
设当有(i)个方块时，记红绿都为偶数的方案数为(a_i),一个位偶数的方案数为(b_i),都不是的方案数为(c_i)。
容易得到：

[egin{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} ]

那么构造：

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

方法二：
学过生成函数的同学都知道这是指数生成函数的入门题。
容易得到：

[egin{aligned}g^{(e)}&=(1+dfrac{x^2}{2!}+dfrac{x^4}{4!}+......)^2(1+dfrac{x}{1!}+dfrac{x^2}{2!}+......)^2cr&=(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} ]

容易得到递推式：

[egin{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)，构造：

[egin{pmatrix}w_{i-1}\w_{i-2}\w_{i-3}\2i\2end{pmatrix}×egin{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&1end{pmatrix} ]

我们再运用拓展欧拉定理，来优化系数：不会的戳这里
由于(varphi(1e9+7)=1e9+6),那我们把系数(mod 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
",ans%MODD);
	return 0;
}

这种模型的原理：
利用已知的数列元素，把后面每一项表示出来，这是通常是这个元素幂的积的形式，通常使用拓展欧拉定理优化幂次。
核心思想是：同底数幂相乘，底数不变，指数相加。

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

【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/tlx-blog/p/12676643.html