[原创]HDU 5895&&2016 ACM/ICPC Asia Regional Shenyang Online1004 Mathematician QSC [矩阵加速+欧拉降幂]【数论】

# include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxn = 505;
const double  Pi =  acos(-1);
# define pb push_back
# define lalal puts("****");
const int M = 4;
int MOD ;
struct Matrix
{
    LL m[M][M];
    void clearO()
    {
        for(int i=0; i<M; i++) //初始化矩阵
            for(int j=0; j<M; j++)
                m[i][j]= 0;
    }
    void clearE()
    {
        for(int i=0; i<M; i++) //初始化矩阵
            for(int j=0; j<M; j++)
                m[i][j]= (i==j);
    }
    void display()
    {
        for(int i=0; i<M; i++)
            {
                for(int j=0; j<M; j++)
                printf("%d ",m[i][j]);
                puts("");
            }
    }
};

Matrix operator * (Matrix a,Matrix b)
{
    Matrix c;
    c.clearO();

    for(int k=0; k<M; k++)
        for(int i=0; i<M; i++) //实现矩阵乘法
        {
            if(a.m[i][k] <= 0)  continue;
            for(int j=0; j<M; j++)
            {
                if(b.m[k][j] <= 0)    continue;
                c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]+MOD)%MOD;
            }
        }
    return c;
}

Matrix operator ^ (Matrix a,LL b)
{
    Matrix c;
    c.clearE();
    while(b)
    {
        if(b&1) c= c * a ;
        b >>= 1;
        a = a * a ;
    }
    return c;
}

int Is_or[101001];
int prime[13000],kpri;
void Prime()
{
    int n=100001;
    kpri=0;
    memset(Is_or,1,sizeof(Is_or));
    Is_or[0]=Is_or[1]=0;
    for(int i=2; i<n; i++)
    {
        if(Is_or[i])
        {
            prime[kpri++]=i;
            for(int j=i+i; j<n; j+=i)
            {
                Is_or[j]=0;
            }
        }
    }
    //prime[kpri]=10007;
    //printf("%d\n",kpri);
    return ;
}

LL Phi(LL n)
{
    LL rea=n;
    for(int i=0; prime[i]*prime[i]<=n; i++)
    {
        if(n%prime[i]==0)
        {
            rea=rea-rea/prime[i];
            while(n%prime[i]==0)
                n/=prime[i];
        }
    }
    if(n>1)  rea=rea-rea/n;
    return rea;
}

LL qmod(LL a,LL b)
{
    LL res= 1;
    while(b)
    {
        if(b&1) res=(res*a)%MOD;
        b>>=1;
        a=(a*a)%MOD;
    }
    return res;
}

int main()
{
    Prime();
    int _;
    while(~scanf("%d",&_))
    {
        while(_--)
        {
            LL n,y,x,s;
            scanf("%I64d%I64d%I64d%I64d",&n,&y,&x,&s);

            MOD = Phi(s+1);

            Matrix a,b;
            a.clearO(),b.clearO();
            a.m[0][1]=1;

            b.m[0][1]=1;
            b.m[1][0]=1,b.m[1][1]=4,b.m[1][2]=2,b.m[1][3]=1;
            b.m[2][1]=4,b.m[2][2]=1;
            b.m[3][3]=1;

            b=b^(n*y);
            a=a*b;

            LL zhi = a.m[0][3]%MOD+MOD;

            MOD=s+1;
            printf("%I64d\n",qmod(x%MOD,zhi));
        }
    }
    return 0;
}