好题!
/*
gi=c^i * fi
gi=gi-1 * gi-2 * gi-3
把 g1,g2,g3 质因数分解
g1=p1^e11 * p2^e12 * p3^e13 ... pk^e1k
g2=p1^e21 * p2^e22 * p3^e23 ... pk^e2k
g3=p1^e31 * p2^e32 * p3^e33 ... pk^e3k
然后构造初始矩阵
e11 e12 e13 ... e1k
e21 e22 e23 ... e2k
e31 e32 e33 ... e3k
这个 3*m 矩阵前面乘以系数矩阵 3*3
0 1 0
1 0 0
1 1 1
矩阵快速幂得到最后矩阵
en1 en2 en3 ... enk
...
...
即每个质因子最后的指数
然后计算出 gn=c^n * fn
fn = gn*inv(c^n)
注意!!! a^x = a^(x % phi(p)) (mod p)
所以矩阵里乘法的模数时 mod-1 !
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define maxn 2000005
#define mod 1000000007
ll n,c,f1,f2,f3,A[50][50],M[50][50];
set<ll>s;
ll p[50],E[50][50],k;
void divide(ll x){
ll tmp=x;
for(ll i=2;i*i<=tmp;i++)
if(tmp%i==0){
s.insert(i);
while(tmp%i==0)tmp/=i;
}
if(tmp>1)s.insert(tmp);
}
void calc(ll x,int id){
for(ll i=1;i<=k;i++)
while(x%p[i]==0)
E[id][i]++,x/=p[i];
}
void mul1(ll A[50][50],ll B[50][50]){
ll res[50][50]={};
for(int i=1;i<=3;i++)
for(int j=1;j<=3;j++)
for(int k=1;k<=3;k++)
res[i][j]=(res[i][j]+A[i][k]*B[k][j]%(mod-1))%(mod-1);
memcpy(A,res,sizeof res);
}
void mul2(ll A[50][50],ll B[50][50]){
ll res[50][50]={};
for(int i=1;i<=3;i++)
for(int j=1;j<=k;j++)
for(int l=1;l<=3;l++)
res[i][j]=(res[i][j]+A[i][l]*B[l][j]%(mod-1))%(mod-1);
memcpy(B,res,sizeof res);
}
ll Pow(ll a,ll b){
ll res=1;
while(b){
if(b%2)
res=res*a%mod;
b>>=1;a=a*a%mod;
}
return res;
}
int main(){
cin>>n>>f1>>f2>>f3>>c;
divide(f1);divide(f2);divide(f3);divide(c);
for(auto it:s)
p[++k]=it;
calc(f1*c,1);
calc(f2*c,2);calc(c,2);
calc(f3*c,3);calc(c,3);calc(c,3);
A[1][1]=0;A[1][2]=1;A[1][3]=0;
A[2][1]=0;A[2][2]=0;A[2][3]=1;
A[3][1]=1;A[3][2]=1;A[3][3]=1;
ll b=n-1,res[50][50]={};
for(int i=1;i<=3;i++)
res[i][i]=1;
while(b){
if(b%2)
mul1(res,A);
b>>=1;mul1(A,A);
}
mul2(res,E);
ll ans=1;//ans=c^n * fn
for(int i=1;i<=k;i++)
ans=ans*Pow(p[i],E[1][i])%mod;
ll tmp=Pow(c,n);
tmp=Pow(tmp,mod-2);// 求逆元
ans=ans*tmp%mod;
cout<<ans<<'\n';
}
*/
来源: http://www.bubuko.com/infodetail-3090620.html