HDU3221Brute-force Algorithm(矩陣快速冪&&指數降冪)

題目連結：http://acm.hdu.edu.cn/showproblem.php?pid=3221

Professor Brute is not good at algorithm design. Once he was asked to solve a path finding problem. He worked on it for several days and finally came up with the following algorithm:

Any fool but Brute knows that the function “funny” will be called too many times. Brute wants to investigate the number of times the function will be called, but he is too lazy to do it.

Now your task is to calculate how many times the function “funny” will be called, for the given a, b and n. Because the answer may be too large, you should output the answer module by P.

從第三項開始 a的指數為斐波那契數列

由於比較大 我們根據公式a^b%p=a^(b%phi(p)+phi(p))%p   b>=phi(p)對指數進行降冪

然後再快速冪取模即可得到結果

對p=1的情況進行特判；

#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;
typedef long long LL;
const int N = 2;
struct matrax
{
   LL m[N][N];
};
matrax E={
 1,0,
 0,1
};
matrax A={
 1,1,
 1,0
};
LL a,b,p,n,mm;
matrax multi(matrax a,matrax b)
{
    matrax c;
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            c.m[i][j]=0;
            for(int k=0;k<N;k++)
                c.m[i][j]+=a.m[i][k]*b.m[k][j];
            if(c.m[i][j]>mm)
                c.m[i][j]=c.m[i][j]%mm+mm;
        }
    }
    return c;
}
matrax pow(matrax a,LL n)
{
    matrax ans=E,p=a;
    while(n){
        if(n&1){
            ans=multi(ans,p);
            n--;
        }
        n>>=1;
        p=multi(p,p);
    }
    return ans;
}
LL phi(LL n)
{
    LL rea=n;
    for(int i=2;i*i<=n;i++){
        if(n%i==0){
            rea=rea-rea/i;
            while(n%i==0)
                n/=i;
        }
    }
    if(n>1)
        rea=rea-rea/n;
    return rea;
}
LL quick_mod(LL a,LL b)
{
    LL ans=1;
    while(b){
        if(b&1){
            ans=ans*a%p;
            b--;
        }
        b>>=1;
        a=a*a%p;
    }
    return ans;
}
int main()
{
    int t=1,cas;
    cin>>cas;
    while(cas--){
        cin>>a>>b>>p>>n;
        mm=phi(p);
        printf("Case #%d: ",t++);
        if(n==1){
            printf("%I64d\n",a%p);
            continue;
        }
        if(n==2){
            printf("%I64d\n",b%p);
            continue;
        }
        if(n==3){
            printf("%I64d\n",a*b%p);
            continue;
        }
        if(p==1){
            puts("0");
            continue;
        }
        matrax g=pow(A,n-2);
        LL m1,m2,num1,num2;
        m1=g.m[1][0];
        m2=g.m[0][1]+g.m[1][1];
        if(m2>mm) m2=m2%mm+mm;
        num1=quick_mod(a,m1);
        num2=quick_mod(b,m2);
        printf("%I64d\n",num1*num2%p);
    }
    return 0;
}
#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;
typedef long long LL;
const int N = 2;
struct matrax
{
   LL m[N][N];
};
matrax E={
 1,0,
 0,1
};
matrax A={
 1,1,
 1,0
};
LL a,b,p,n,mm;
matrax multi(matrax a,matrax b)
{
    matrax c;
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            c.m[i][j]=0;
            for(int k=0;k<N;k++)
                c.m[i][j]+=a.m[i][k]*b.m[k][j];
            if(c.m[i][j]>mm)
                c.m[i][j]=c.m[i][j]%mm+mm;
        }
    }
    return c;
}
matrax pow(matrax a,LL n)
{
    matrax ans=E,p=a;
    while(n){
        if(n&1){
            ans=multi(ans,p);
            n--;
        }
        n>>=1;
        p=multi(p,p);
    }
    return ans;
}
LL phi(LL n)
{
    LL rea=n;
    for(int i=2;i*i<=n;i++){
        if(n%i==0){
            rea=rea-rea/i;
            while(n%i==0)
                n/=i;
        }
    }
    if(n>1)
        rea=rea-rea/n;
    return rea;
}
LL quick_mod(LL a,LL b)
{
    LL ans=1;
    while(b){
        if(b&1){
            ans=ans*a%p;
            b--;
        }
        b>>=1;
        a=a*a%p;
    }
    return ans;
}
int main()
{
    int t=1,cas;
    cin>>cas;
    while(cas--){
        cin>>a>>b>>p>>n;
        mm=phi(p);
        printf("Case #%d: ",t++);
        if(n==1){
            printf("%I64d\n",a%p);
            continue;
        }
        if(n==2){
            printf("%I64d\n",b%p);
            continue;
        }
        if(n==3){
            printf("%I64d\n",a*b%p);
            continue;
        }
        if(p==1){
            puts("0");
            continue;
        }
        matrax g=pow(A,n-2);
        LL m1,m2,num1,num2;
        m1=g.m[1][0];
        m2=g.m[0][1]+g.m[1][1];
        if(m2>mm) m2=m2%mm+mm;
        num1=quick_mod(a,m1);
        num2=quick_mod(b,m2);
        printf("%I64d\n",num1*num2%p);
    }
    return 0;
}