HDU 5895 Mathematician QSC

矩阵快速幂，欧拉定理。

$g(n)$递推式：$g(n)=5g(n-1)+5g(n-2)-g(n-3)$，可以构造矩阵快速求递$n$项，指数很大，可以利用欧拉定理降幂。

#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<iostream>
using namespace std;
typedef long long LL;
const double pi=acos(-1.0),eps=1e-6;
void File()
{
    freopen("D:\\in.txt","r",stdin);
    freopen("D:\\out.txt","w",stdout);
}
template <class T>
inline void read(T &x)
{
    char c=getchar(); x=0;
    while(!isdigit(c)) c=getchar();
    while(isdigit(c)) {x=x*10+c-'0'; c=getchar();}
}

LL MOD;

struct Matrix
{
    long long A[4][4];
    int R, C;
    Matrix operator*(Matrix b);
};

Matrix X, Y, Z;

Matrix Matrix::operator*(Matrix b)
{
    Matrix c;
    memset(c.A, 0, sizeof(c.A));
    int i, j, k;
    for (i = 1; i <= R; i++)
        for (j = 1; j <= b.C; j++)
            for (k = 1; k <= C; k++)
                c.A[i][j] = (c.A[i][j] + (A[i][k] * b.A[k][j]) % MOD) % MOD;
    c.R = R; c.C = b.C;
    return c;
}

void init()
{
    memset(X.A, 0, sizeof X.A);
    memset(Y.A, 0, sizeof Y.A);
    memset(Z.A, 0, sizeof Z.A);

    Y.R = 3; Y.C = 3;
    for (int i = 1; i <= 3; i++) Y.A[i][i] = 1;

    X.R = 3; X.C = 3;
    X.A[1][3]=-1;
    X.A[2][1]=1; X.A[2][3]=5;
    X.A[3][2]=1; X.A[3][3]=5;

    Z.R = 1; Z.C = 3;
    Z.A[1][1]=0; Z.A[1][2]=1; Z.A[1][3]=5;
}

LL work(LL tt)
{
    if(tt==0) return 0;
    if(tt==1) return 1;
    if(tt==2) return 5;

    tt=tt-2; init();
    while (tt)
    {
        if (tt % 2 == 1) Y = Y*X;
        tt = tt >> 1;
        X = X*X;
    }
    Z = Z*Y;
    return Z.A[1][3];
}

int T;
LL n,x,y,s;

LL phi(LL x)
{
     LL res=x,a=x;
     for(int i=2;i*i<=a;i++)
     {
         if(a%i==0)
         {
             res=res/i*(i-1);
             while(a%i==0) a/=i;
         }
     }
     if(a>1) res=res/a*(a-1);
     return res;
}

LL pow(LL a,LL b,LL mod)
{
    LL c=1;
    while(b!=0)
    {
        if(b%2==1) c=(c*a)%mod,b--;
        else a=(a*a)%mod,b=b/2;
    }
    return c;
}

int main()
{
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%lld%lld%lld",&n,&y,&x,&s);
        MOD=phi(s+1);
        LL gn=work(n*y)+MOD;
        LL ans=pow(x,gn,s+1);
        printf("%lld\n",ans);
    }
    return 0;
}