Brute-force Algorithm [HDU 3221]

http://acm.hdu.edu.cn/showproblem.php?pid=3221

对于指数取欧拉函数需要注意一下,a^x%P=a^(x%phi(P)+phi(P))%P

//A^x = A^(x % Phi(C) + Phi(C)) (mod C) 的应用
const int MM = 1111111;
#define debug puts("wrong")
typedef __int64 int64;
const double lep=1e-10;
int64 a,b,P,n,mod;

const int maxn = 3; //矩阵大小
struct Matrix {
    int64 r, c; //矩阵大小
    int64 a[maxn][maxn];
    //初始化
    Matrix(int64 r,int64 c):r(r), c(c) {}
    void reset() {memset(a,0,sizeof(a));}
    //重写(),方便存取
    int64& operator() (int64 i,int64 j) {return a[i][j];}
    //O(N^3) 矩阵乘法
    Matrix operator*(const Matrix&B) {
        Matrix M(r,B.c); M.reset();
        for(int i=0;i<r;i++) {
            for(int j=0;j<B.c;j++) {
                for(int k=0;k<c;k++) {
                    M.a[i][j]=M.a[i][j]+a[i][k]*B.a[k][j];
                    if(M.a[i][j]>mod) M.a[i][j]=M.a[i][j]%mod+mod; //以防出现可以整除的情况
                }
            }
        }
        return M;
    }
    Matrix operator^(int64 n) {
        Matrix M(r,c); M.reset();
        for(int i=0;i<r;i++) M(i,i)=1;
        Matrix A=*this;
        for(;n;n>>=1,A=A*A) if(n&1) M=A*M;
        return M;
    }
};

int64 pow(int64 x,int64 y) {
    int64 res=1;
    while(y) {
        if(y&1) {
            res=res*x;
            if(res>=P) res%=P;
        }
        x=x*x; y>>=1;
        if(x>=P) x%=P;
    }
    return res;
}
int64 euler(int64 x) {
    int64 i,res=x;
    for(i=2;(i*i)<=x;i++) {
        if(x%i==0) {
            res=res/i*(i-1);
            while(x%i==0) x/=i;
        }
    }
    if(x>1) res=res/x*(x-1);
    return res;
}
void get_data() {
    int i,j,k;
    scanf("%I64d%I64d%I64d%I64d",&a,&b,&P,&n);
    a%=P; b%=P;
}

void solve() {
    int64 i,j,k,res=1,x,y;
    
    printcase();
    if(n==1) {printf("%I64d\n",a);return;}
    if(n==2) {printf("%I64d\n",b);return;}

    Matrix tmp(2,2),tt(2,1),ans(2,2);
    mod=euler(P);
    tmp(0,0)=tmp(0,1)=tmp(1,0)=1, tmp(1,1)=0; 
    tmp=tmp^(n-2); 
    
    tt(0,0)=1, tt(1,0)=0;
    ans=tmp*tt; y=ans(0,0);

    tt(0,0)=0, tt(1,0)=1;
    ans=tmp*tt; x=ans(0,0);
    
    res=(res*pow(a,x))%P;
    res=(res*pow(b,y))%P;
    printf("%I64d\n",res);
}
int main() {
    int ca; scanf("%d",&ca);
    while(ca--) get_data(),solve();
    return 0;
}