hdu 3221 Brute-force Algorithm

http://acm.hdu.edu.cn/showproblem.php?pid=3221
ans=a^f 1[n]*b^f 2[n]を要約することができます.
ここで、f 1[1]=1,f 1[2]=0,f 1[n]=f 1[n-1]+f 1[n-2]
            f 2[1]=0,f 2[2]=1,f 2[n]=f 2[n-1]+f 2[n-2]
数式に従ってもいいです  A^x  % C=A^(x%phi(C)+phi(C))%C   x満足x>=phi[C]
それから、マトリックスで素早くべき乗してf 1[n]%phi[p]を求めます.   和f 2[n]%phi[p]   nが大きいので、直接に頼むのは無理です.

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>

using namespace std;

typedef unsigned long long ll;
ll mod;
ll get_eulr(ll a){
    ll ret=a;
    for(int i=2;i*i<=a;i++)
      if(a%i==0){
         ret=ret/i*(i-1);
         while(a%i==0) a/=i;
      }
   if(a!=1) ret=ret/a*(a-1);
   return ret;
}
ll quick_mod(ll a,ll n,ll mo){
     if(!a) return 0;
     ll ans=1;
     while(n){
         if(n&1) ans=ans*a%mo;
         n>>=1;
         a=a*a%mo;
     }
     return ans;
}
void mul(ll a[][2],ll b[][2],ll c[][2]){
    for(int i=0;i<2;i++)
    for(int j=0;j<2;j++){
        c[i][j]=0;
        for(int r=0;r<2;r++){
           c[i][j]=(c[i][j]+a[i][r]*b[r][j])%mod;
        }
    }
}
ll ff(int id,ll n){
   ll a[2][2],b[2][2],c[2][2];
   a[0][0]=0;
   a[0][1]=a[1][0]=a[1][1]=1;
   c[0][0]=c[1][1]=1;
   c[0][1]=c[1][0]=0;
   while(n){
       if(n&1){
           mul(c,a,b);
           memcpy(c,b,sizeof(b));
       }
       n>>=1;
       mul(a,a,b);
       memcpy(a,b,sizeof(b));
   }
   if(!id) return c[0][0];
   else return c[0][1];
}
ll fib[110];
int main()
{
    int cas=1,ca;
    scanf("%d",&ca);
    ll p,n,d1,d2,a,b,ans1,ans2;
    while(ca--)
    {
        scanf("%I64u%I64u%I64u%I64u",&a,&b,&p,&n);
        ll phi=get_eulr(p);
        mod=phi;
        fib[1]=1,fib[2]=0;
        if(n==1) d1=1;
        else if(n==2) d1=0;
        else{
         int i;
         for(i=3;i<=n;i++){
            fib[i]=fib[i-1]+fib[i-2];
            if(fib[i]>=phi) break;
         }
         if(i<=n) d1=ff(0,n-1)+phi;
         else d1=fib[n];
        }

        fib[1]=0,fib[2]=1;
        if(n==1) d2=0;
        else if(n==2) d2=1;
        else{
          int i;
          for(i=3;i<=n;i++){
             fib[i]=fib[i-1]+fib[i-2];
             if(fib[i]>=phi) break;
          }
          if(i<=n) d2=ff(1,n-1)+phi;
          else d2=fib[n];
        }
        ans1=quick_mod(a,d1,p);
        ans2=quick_mod(b,d2,p);
        printf("Case #%d: %I64u
",cas++,ans1*ans2%p);
    }
    return 0;
}