NYOJ 1000

整体的思路是“快速幂取模算法”和“矩阵快速幂算法”的结合

F[0] = a¹b⁰

F[1] = a⁰b¹

F[2] = a¹b¹

F[3] = a¹b²

F[4] = a²b³

F[5] = a³b⁵

F[6] = a⁵b⁸

...

观察可知F[n] = a^f(n-2)·b^f(n-1)其中f(n-1)、f(n-2)是斐波那契数列，设p=1000000007(经检验可知p是质数,这引导我们中间计算矩阵幂时可用费马小定理)

故问题可转化为求  a^f(n-2)·b^f(n-1)%p，即 (a^f(n-2)%p)·(b^f(n-1)%p)%p.  接下来用快速幂取模算法解决，详见 http://www.cnblogs.com/liugl7/p/4831030.html

在使用快速幂取模算法之前需要解决的问题是快速计算斐波那契数列的值的问题，传统的迭代算法并不能满足快速计算的需求，所以这里使用矩阵快速幂算法。

令：

                    

则易得  M_n  = M₀Fⁿ   (n>=1)

所以只需要计算出Fⁿ然后乘上M₀即可得到M_n, 然后就可以得到f(n-1),f(n-2)代入  (a^f(n-2)%p)·(b^f(n-1)%p)%p 的计算。

代码如下：

#include <iostream> 
#include <cstdio>  
using namespace std;  
typedef long long LL; 
#define MOD 1000000007  

struct Matrax {     //定义2*2的矩阵 
    LL mat[2][2]; 
}U,F;  

//2*2的两矩阵的乘法 
Matrax multi(Matrax a,Matrax b) {
    Matrax ans;     
    int i,j,k;     
    for(i=0;i<2;i++) {         
        for(j=0;j<2;j++) {             
            ans.mat[i][j]=0;             
            for(k=0;k<2;k++)     
                ans.mat[i][j]=(ans.mat[i][j]+(a.mat[i][k]%(MOD-1))*(b.mat[k][j]%(MOD-1)))%(MOD-1); 
         }      
    }         
    return ans;  
}  

//实现矩阵p的n次方   
Matrax my_pow_matrax(int n) {     
    Matrax p=F,ans=U;     
    while(n) {         
        if(n&1) ans=multi(ans,p);
                 
        n>>=1;         
        p=multi(p,p);     
    }     
    return ans; 
}  

//初始化矩阵U=[1,0;0,1]，即二阶单位矩阵；  F=[1,1;1,0] 
void Init() {     
    U.mat[0][0]=U.mat[1][1]=1;     
    U.mat[0][1]=U.mat[1][0]=0;          
    F.mat[0][0]=F.mat[0][1]=F.mat[1][0]=1;     
    F.mat[1][1]=0; 
}  

//快速幂取模算法,可用来得到  a^f(n-2)%p  和  b^f(n-1)%p   
LL my_pow_ll(LL a,LL n) {     
    LL ans=1ll,p=a;     
    while(n)     {         
        if(n&1)             
            ans=((ans%MOD)*(p%MOD))%MOD;         
        n>>=1;         
        p=((p%MOD)*(p%MOD))%MOD;      
    }     
    return ans; 
}  

int main() {     
    Init();     
    LL a,b,n;     
    while(scanf("%lld%lld%lld",&a,&b,&n)!=EOF) { 
        if(n==0) {             
            printf("%lld
",a%MOD);
             continue;         
        }         
        else if(n==1) { 
            printf("%lld
",b%MOD);
             continue;         
        }
         else if(n==2) { 
            printf("%lld
",a*b%MOD);
             continue;         
        }
         else         {             
             Matrax temp=my_pow_matrax(n-1);
             LL f1=temp.mat[0][0]*1+temp.mat[1][0]*0; // fib(n-1)
             LL f2=temp.mat[0][1]*1+temp.mat[1][1]*0;  // fib(n-2)
             LL ans=(my_pow_ll(a,f2)%MOD)*(my_pow_ll(b,f1)%MOD)%MOD;
             printf("%lld
",ans);         
        }     
    }     

    return 0;
}