Codeforces1182E Product Oriented Recurrence(递推+矩乘快速幂）

题目链接
这道题是道好题
给定柿子$f_x=c^{2x-6}\cdot f_{x-1}\cdot f_{x-2}\cdot f_{x-3}c,f_{x-1},f_{x-2},f_{x-3}$已知
求$f{}_n,n\le 10^9$
这题一看就是矩阵快速幂…考虑怎么构造
传统的矩阵乘法只资瓷前n项相加的形式。考虑变形，发现我们可以变成
$f_x\cdot c^x=c^{x-1}\cdot f_{x-1}\cdot c^{x-2}\cdot f_{x-2}\cdot c^{x-3}\cdot f_{x-3}$
设$g_x=f_x\cdot c^{2x}$则这道题变为简单递推
发现$f{}_n\le 10^9$这也就是说这个数质因子个数不超过30个，我们分别求出各个数的质因子，那么数的乘法就是质因子的相加
最后求个逆元回推就好了
这道题用到的两个思想都是非常常见的，下次做题要优先考虑
CF评论区里有个神人，直接把$c$放进矩阵了…orz
Code

#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
struct mat
{
	long long a[5][5];
	mat operator *(const mat z)
	{
		mat f;
		memset(f.a,0,sizeof(f.a));
		for(int i=1;i<=3;i++)
		for(int j=1;j<=3;j++)
		for(int k=1;k<=3;k++)
		f.a[i][j]=(f.a[i][j]+a[k][j]*z.a[i][k]%(mod-1))%(mod-1);
		return f;
	}
	friend ostream &operator <<(ostream &output,const mat &qq)
	{
		for(int i=1;i<=3;i++,output<<endl)
		for(int j=1;j<=3;j++)
		output<< qq.a[i][j]<<" ";
	}
}e,prime[500];
int p[500];
int o;
long long n,f1,f2,f3,c;
map<int,int>mp;
int num[4][55];
void work(int nm,long long x)
{
	mat dd=e;
	while(x)
	{
	//	cout<<prime[nm]<<endl<<dd<<endl;
		if(x&1)prime[nm]=prime[nm]*dd;
		//cout<<prime[nm]<<endl;	
	//	cout<<"fdsasaesf";
		dd=dd*dd;
		x>>=1;
	}
}
long long ksm(long long x,long long a)
{
	long long re=1;
	while(a)
	{
		if(a&1)re=re*x%mod;
		x=x*x%mod;
		a>>=1;
	}
	return re;
}
long long exgcd(long long a,long long b,long long &x,long long &y)
{
	if(b==0)
	{
		x=1;
		y=0;
		return a;
	}
	int t=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return t;
}
int main()
{
	e.a[1][2]=e.a[2][3]=e.a[3][1]=e.a[3][2]=e.a[3][3]=1;
	scanf("%lld%lld%lld%lld%lld",&n,&f1,&f2,&f3,&c);
	f1=f1*c%mod;
	f2=f2*c%mod*c%mod;
	f3=f3*c%mod*c%mod*c%mod;
	long long k1=f1,k2=f2,k3=f3;
	for(int i=2;i<=sqrt(f1);i++)
	{
		if(k1%i==0)
		{
			p[++o]=i;
			while(k1%i==0)
			{
				k1/=i;
			}
		}
		if(k1==1)break;
	}
	if(k1!=1)p[++o]=k1;
	for(int i=2;i<=sqrt(f2);i++)
	{
		if(k2%i==0)
		{
			p[++o]=i;
			while(k2%i==0)
			{
				k2/=i;
			}
		}
		if(k2==1)break;
	}
	if(k2!=1)p[++o]=k2;
	for(int i=2;i<=sqrt(f3);i++)
	{
		if(k3%i==0)
		{
			p[++o]=i;
			while(k3%i==0)
			{
				k3/=i;
			}
		}
		if(k3==1)break;
	}
	if(k3!=1)p[++o]=k3;
	sort(p+1,p+1+o);
	o=unique(p+1,p+1+o)-p-1;
	for(int i=1;i<=o;i++)
	{
		mp[p[i]]=i;
	}
	k1=f1,k2=f2,k3=f3;
	for(int i=2;i<=sqrt(f1);i++)
	{
		if(k1%i==0)
		{
			while(k1%i==0)
			{
				k1/=i;num[1][mp[i]]++;
			}
		}
		if(k1==1)break;
	}
	if(k1!=1)num[1][mp[k1]]++;
	for(int i=2;i<=sqrt(f2);i++)
	{
		if(k2%i==0)
		{
			while(k2%i==0)
			{
				k2/=i;num[2][mp[i]]++;
			}
		}
		if(k2==1)break;
	}
	if(k2!=1)num[2][mp[k2]]++;
	for(int i=2;i<=sqrt(f3);i++)
	{
		if(k3%i==0)
		{
			while(k3%i==0)
			{
			//	cout<<i<<endl;
				k3/=i;num[3][mp[i]]++;
			}
		}
		if(k3==1)break;
	}
	if(k3!=1)num[3][mp[k3]]++;
//	cout<<num[3][2]<<"---";
	long long ans=1ll;
	for(int i=1;i<=o;i++)
	{
	//	cout<<p[i]<<"fds"<<endl;
		prime[i].a[1][1]=num[1][i];
		prime[i].a[2][1]=num[2][i];
		prime[i].a[3][1]=num[3][i];
	//	cout<<prime[i].a[1][1]<<" "<<prime[i].a[2][1]<<" "<<prime[i].a[3][1]<<endl;
		work(i,n-3);
	//	cout<<prime[i].a[1][1]<<" "<<prime[i].a[2][1]<<" "<<prime[i].a[3][1]<<"fds"<<endl;
		ans=ans*ksm(p[i],prime[i].a[3][1])%mod;
	//	cout<<ans<<endl;
	}
	
	long long modp,y;
	c=ksm(c,n);
	exgcd(c,mod,modp,y);//cout<<"??";
	modp=(modp%mod+mod)%mod;
	printf("%lld",ans*modp%mod);
	return 0;
}