hdu 3221 矩阵幂乘 + 超大a*b%c

for(int n1=0;n1<A.size();n1++){
for(int m2=0;m2<B[0].size();m2++){
for(int m1=0;m1<A[0].size();m1++){
C[n1][m2]=(C[n1][m2]+A[n1][m1]*B[m1][m2]%M)%M;
}
}
}

#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<vector>
#include<cmath>
#include<iostream>
using namespace std;

const long long maxn=100000*5;
typedef long long ll;

typedef vector<ll> vec;
typedef vector<vec> mat;
ll M;
mat mul(mat &A,mat &B){
mat C(A.size(),vec(B[0].size()));

for(int n1=0;n1<A.size();n1++){
for(int m2=0;m2<B[0].size();m2++){
for(int m1=0;m1<A[0].size();m1++){
C[n1][m2]=(C[n1][m2]+A[n1][m1]*B[m1][m2]);
if(C[n1][m2]>M){
C[n1][m2]=(C[n1][m2]%M)+M;
}
}
}
}
return C;
}
mat pow(mat A,ll n){
mat B(A.size(),vec(A.size()));
for(int i=0;i<A.size();i++){
B[i][i]=1;
}
while(n>0ll){
if(n&1ll) B=mul(B,A);
A=mul(A,A);
n>>=1ll;
}
return B;
}
typedef pair<ll,ll> P;
P cal(ll n){
P p;
if(n==1){
p.first=1;p.second=0;

}else if(n==2){
p.first=0;p.second=1;
}else if(n==3){
p.first=1;p.second=1;
}else{

mat A(2,vec(2));
mat B(2,vec(2));
{
A[0][0]=1; A[0][1]=1;
A[1][0]=1; A[1][1]=0;
}

{
B[0][0]=1; B[0][1]=1;
B[1][0]=1; B[1][1]=0;
}
A=pow(A,n-3ll);
B=pow(B,n-2ll);
ll x= A[0][0]*1ll+A[0][1]*0ll;
ll y= B[0][0]*1ll+B[0][1]*0ll;
p.first=x;p.second=y;
}
return p;
}

ll powm(ll x,ll n,ll pmod){
ll res=1;
while(n>0ll){
if(n&1ll)
res=res*x%pmod;

x=x*x%pmod;
n>>=1ll;
}
return res;
}

ll getol(ll n){

ll ans=n;
for(ll i=2;i*i<=n;i++){
if(n%i==0ll){
ans-=ans/i;
while(n%i==0ll){
n/=i;
}
}
}
if(n>1ll) ans-=ans/n;
return ans;
}

int main()
{
// freopen("D://in.txt","r",stdin);
int cas;
cin>>cas;
ll a,b,n,pmod;
for(int ik=1;ik<=cas;ik++){
cin>>a>>b>>pmod>>n;

M=getol(pmod);
P p=cal(n);
ll res;

res=powm(a,p.first,pmod)%pmod;
res=(res*(powm(b,p.second,pmod)%pmod))%pmod;

//  printf("%I64d %I64d %I64d\n",M ,p.first,p.second);
printf("Case #%d: %I64d\n",ik,res);
}

return 0;
}