HDU 5895 Mathematician QSC

矩阵快速幂，欧拉函数

Hint: \(g(n)=\frac{f(n)f(n+1)}{2}\)

#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,b) for(int i=a;i<=b;++i)
#define ms(arr,a) memset(arr,a,sizeof arr)
#define debug(x) cout<<"< "#x" = "<<x<<" >"<<endl

long long n,y,x,s;
long long F;

struct matrix
{
    long long a,b,c,d;
    matrix(){}
    matrix(long long m,long long n,long long p,long long q):a(m),b(n),c(p),d(q){}
    matrix(long long x):a(x),b(0),c(0),d(x){}
};
matrix solve(matrix A,matrix B)
{
    matrix ret;
    ret.a=(A.a*B.a+A.b*B.c)%F;
    ret.b=(A.a*B.b+A.b*B.d)%F;
    ret.c=(A.c*B.a+A.d*B.c)%F;
    ret.d=(A.c*B.b+A.d*B.d)%F;
    return ret;
}
long long solve(long long x,long long y)
{
    return x*y%(s+1);
}
template<typename T>
T quick_pow(T a,long long n)
{
    T ret(1);
    while(n)
    {
        if(n&1)ret=solve(ret,a);
        a=solve(a,a);
        n>>=1;
    }
    return ret;
}
long long phi(long long n)
{
    long long ret=n;
    for(long long i=2;i<(long long)sqrt(n)+1;++i)
    {
        if(n%i==0)ret=ret/i*(i-1);
        while(n%i==0)n/=i;
    }
    if(n>1)ret=ret/n*(n-1);
    return ret;
}
int main()
{
    int T;scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%lld%lld%lld",&n,&y,&x,&s);
        F=2*phi(s+1);
        long long N=n*y;
        matrix M(2,1,1,0);
        M=quick_pow(M,N);
        printf("%lld\n",quick_pow(x,M.a*M.c%F/2));
    }
}