CF #566 Div. 2 E Product Oriented Recurrence//矩阵快速幂+费马小定理

题意：

 思路：

最后f(n)=f(3)^[0][0] * f(2)^[0][1] * f(1)^[0][2] * c^g(n)

费马小定理：如果p是一个质数，而整数a不是p的倍数，则有a^（p-1）≡1（mod p）。

然后这里求的都是系数，系数取模需要mod1e9+6。（用到费马小定理，a^(1e9+7-1) ≡ a^0 ）

c的系数是一眼就看出来的 -。-    复杂的想到了，这个把我绕进去了........

#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define FI first
#define SE second
const LL mod =1e9+6;
const int MX = 1e6+5;
LL n,f1,f2,f3,c;
struct mat{
    LL m[5][5];
};
mat matMul(mat aa,mat bb,int sz){
    mat re;memset(re.m,0,sizeof(re.m));
    for(int i=0;i<sz;i++)
        for(int j=0;j<sz;j++)
            for(int k=0;k<sz;k++)
                re.m[i][j]=(re.m[i][j]+aa.m[i][k]*bb.m[k][j])%mod;
    return re;
}
mat matQuickMod(mat jin,LL x,int sz){
    mat re;memset(re.m,0,sizeof(re.m));
    for(int i=0;i<sz;i++) re.m[i][i]=1;
    while(x){
        if(x&1) re=matMul(re,jin,sz);
        x=x>>1;jin=matMul(jin,jin,sz);
    }
    return re;
}
void jinMatInit(mat &x){//乘的转移矩阵
    memset(x.m,0,sizeof(x.m));
    x.m[0][0]=1;x.m[0][1]=1;x.m[0][2]=0;x.m[0][3]=0;x.m[0][4]=0;
    x.m[1][0]=1;x.m[1][1]=0;x.m[1][2]=1;x.m[1][3]=0;x.m[1][4]=0;
    x.m[2][0]=1;x.m[2][1]=0;x.m[2][2]=0;x.m[2][3]=0;x.m[2][4]=0;
    x.m[3][0]=2;x.m[3][1]=0;x.m[3][2]=0;x.m[3][3]=1;x.m[3][4]=0;
    x.m[4][0]=-6;x.m[4][1]=0;x.m[4][2]=0;x.m[4][3]=1;x.m[4][4]=1;
}
void secondJinMatInit(mat &x){
    memset(x.m,0,sizeof(x.m));
    x.m[0][0]=1;x.m[0][1]=1;x.m[0][2]=1;
    x.m[1][0]=1;x.m[1][1]=0;x.m[1][2]=0;
    x.m[2][0]=0;x.m[2][1]=1;x.m[2][2]=0;
}
LL quickPow(LL x,LL y){
    LL re=1;LL jin=x;while(y){
        if(y%2) re=re*jin%(mod+1);
        y=y/2;jin=jin*jin%(mod+1);
    }
    return re;
}
int main(){
    cin>>n>>f1>>f2>>f3>>c;
    mat G;jinMatInit(G);
    mat B=matQuickMod(G,n-3,5);
    mat x;memset(x.m,0,sizeof(x.m));
    x.m[0][0]=0;x.m[0][1]=0;x.m[0][2]=0;x.m[0][3]=4;x.m[0][4]=1;
    x=matMul(x,B,5);
    LL ans=quickPow(c,x.m[0][0]);
    secondJinMatInit(G);
    B=matQuickMod(G,n-4,3);
    memset(x.m,0,sizeof(x.m));
    x.m[0][0]=1;x.m[0][1]=1;x.m[0][2]=1;
    x=matMul(x,B,3);
    ans=ans*quickPow(f3,x.m[0][0])%(mod+1);
    ans=ans*quickPow(f2,x.m[0][1])%(mod+1);
    ans=ans*quickPow(f1,x.m[0][2])%(mod+1);
    cout<<ans;
    return 0;
}