好题！
/* 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<<' '; }
*/