Codeforces Round #566 (Div. 2)-E. Product Oriented Recurrence

地址：https://codeforces.com/contest/1182/problem/E

思路：矩阵快速幂

fn=c^(2n−6)⋅f(n−1)⋅f(n−2)⋅f(n−3)=c^sc * f1^s1 * f2^s2 * f3^s3

因此分别算出s1,s2,s3和sc的值即可算出 fn

而由f4开始按照fn=c^(2n−6)⋅f(n−1)⋅f(n−2)⋅f(n−3)=c^sc * f1^s1 * f2^s2 * f3^s3推

       f1   f2    f3    c

f4    1     1     1    2

f5    1     2     2    6

f6    2     3     4    14

f7    4     6     7    30

f8    7    11    13   60

到f7时 对f1(7):f7=f6+f5+f4   f2(7)=f1(7)+f1(6)    f3(7)=f1(8)    而 fc(7)=fc(6)+fc(5)+fc(4)+2*7-6

因此

f1(n)=f1(n-1)+f1(n-2)+f1(n-3)

f2(n)=f1(n)+f1(n-1)

f3(n)=f1(n+1)

fc(n)=fc(n-1)+fc(n-2)+fc(n-3)+2*n-6

那么就可以用矩阵快速幂来求出f1(n),fc(n)，从而得到 s1,s2,s3,sc

注意的是，在求f1(n),fc(n)时由于 其是 f1,c的幂系数，在取模时要用费马小定理（a^(p-1)=1 mod(p)）来取模 即 MOD=1e9+6

在求f1^s1时 MOD=1e9+7

Code:

#include<iostream>
#include<cstring>
using namespace std;
typedef long long LL;

LL MOD=1e9+6;
const int MAX_S=5;
struct Matrix{
	LL a[MAX_S][MAX_S];
	Matrix (){
		memset(a,0,sizeof(a));
	}
	Matrix operator*(const Matrix &A){
		Matrix B=Matrix();
		for(int k=0;k<MAX_S;++k)
			for(int i=0;i<MAX_S;++i)
				for(int j=0;j<MAX_S;++j)
					B.a[i][j]=(B.a[i][j]+a[i][k]*A.a[k][j])%MOD;
		return B;
	}
};
LL n,c,f[5];

Matrix mypowMa(Matrix A,LL n){
	Matrix res=Matrix();
	if(n<=0)	return res;
	for(int i=0;i<MAX_S;++i)
		res.a[i][i]=1;
	while(n){
		if(n&1)	res=res*A;
		A=A*A;
		n>>=1;
	}
	return res;
}

LL mypowLL(LL a,LL n){
	LL res=1;
	while(n){
		if(n&1)	res=res*a%MOD;
		a=a*a%MOD;
		n>>=1;
	}
	return res;
}
int main()
{
	Matrix A=Matrix(),res;
	A.a[0][0]=A.a[0][1]=A.a[0][2]=1;
	A.a[1][0]=A.a[2][1]=1;
	cin>>n>>f[1]>>f[2]>>f[3]>>c;
	LL s1,s2,s3,sc,ni=n-3,d[5]={0,1,1,2,4};
	if(ni<=3){
		s1=d[ni];	s2=s1+d[ni-1];	s3=d[ni+1];
	}else{
		res=mypowMa(A,ni-3);
		s1=(res.a[0][0]*d[3]%MOD+res.a[0][1]*d[2]%MOD+res.a[0][2]*d[1]%MOD)%MOD;
		res=mypowMa(A,ni-4);
		s2=(s1+(res.a[0][0]*d[3]%MOD+res.a[0][1]*d[2]%MOD+res.a[0][2]*d[1]%MOD)%MOD)%MOD;
		res=mypowMa(A,ni-2);
		s3=(res.a[0][0]*d[3]%MOD+res.a[0][1]*d[2]%MOD+res.a[0][2]*d[1]%MOD)%MOD;
	}
	A.a[0][3]=A.a[0][4]=A.a[3][3]=A.a[3][4]=A.a[4][4]=1;
	res=mypowMa(A,n-3);
	sc=res.a[0][4]*2%MOD;
	MOD+=1;
	LL ans=mypowLL(c,sc);
	ans=ans*mypowLL(f[1],s1)%MOD;
	ans=ans*mypowLL(f[2],s2)%MOD;
	ans=ans*mypowLL(f[3],s3)%MOD;
	cout<<ans<<endl;
	
	return 0;
}