欧拉定理+质因子分解+矩阵快速幂——cf1182E

好题！

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

*/