HDU 5895 矩阵快速幂+欧拉降幂公式+指数循环节

题目链接：HDU5895
注意到题目的数据非常之大，并且递归之后求的指数幂非常之大，所以要先知道欧拉降幂公式:
,相当于形成了指数的循环节。
根据递推关系得到：

建立矩阵:

用矩阵快速幂求解出g(n*y)，之后再利用快速幂取模求出最后的结果

由于数据很大，过程中把数据开大，及时取模。

AC代码：

#include <cstdio>#include <cstring>#include<iostream>#include <cmath>using namespace std;typedef long long LL;typedef unsigned long long ULL;const int N=4;LL n,x,y,s,mod;struct mat{    ULL m[N][N];    mat() {}    mat unit(){        for(int i=0;i<N;i++)            for(int j=0;j<=N;j++)                m[i][j]=i==j?1:0;    }};mat t;mat operator * (mat a,mat b){    //cout<<mod<<endl;    mat res;    for(int i=0;i<N;i++)        for(int j=0;j<=N;j++){            res.m[i][j]=0;            for(int k=0;k<N;k++){                res.m[i][j]+=a.m[i][k]*b.m[k][j];                if(res.m[i][j]>mod)                    res.m[i][j]%=mod;            }    }    return res;}mat operator ^ (mat res,LL n){ //矩阵快速幂    res.unit();    while(n>=1){        if(n&1)            res=res*t;        n=n>>1;        t=t*t;    }    return res;}void init(){    t.m[0][0]=4;t.m[0][1]=4;t.m[0][2]=1;t.m[0][3]=0;    t.m[1][0]=2;t.m[1][1]=1;t.m[1][2]=0;t.m[1][3]=0;    t.m[2][0]=1;t.m[2][1]=0;t.m[2][2]=0;t.m[2][3]=0;    t.m[3][0]=4;t.m[3][1]=4;t.m[3][2]=1;t.m[3][3]=1;}LL quickpow(LL a,LL b,LL c){   //快速幂取模    LL ans=1;    while(b>0){        if(b&1)            ans=ans*a%c;        b=b>>1;        a=(a%c)*(a%c)%c;    }    return ans;}  LL GetEuler(int n)  //欧拉函数{     LL i;     LL res = n,a = n;     for(i = 2;i*i <= a; ++i){         if(a%i == 0){             res -= res/i;              while(a%i == 0)                  a/=i;         }     }     if(a > 1)          res -= res/a;     return res;}int main(){    //freopen("input.txt","r",stdin);    int T;    scanf("%d",&T);    while(T--){        scanf("%lld%lld%lld%lld",&n,&y,&x,&s);//输入顺序搞成nxys，调了两个小时才发现，好气哦！        mat A;        memset(A.m,0,sizeof(A.m));        init();        mod=GetEuler(s+1);        //cout<<mod<<endl;        n=n*y;        if(n==0){            printf("1\n");            continue;        }        else{            LL ans=0,res;            //cout<<n<<endl;            A=t^(n-1);            /*for(int i=0;i<N;i++){                for(int j=0;j<N;j++)                    printf("%lld ",A.m[i][j]);            cout<<endl;            }*/            ans=(A.m[3][0]%mod+A.m[3][3]%mod)%mod;            ans+=mod;            res=quickpow(x,ans,s+1);            printf("%lld\n",res);        }    }    return 0;}