矩阵快速幂+快速幂+欧拉降幂+费马小定理

hdu 4549

M斐波那契数列
Time Limit: 3000/1000 MS (Java/Others) Memory Limit: 65535/32768 K (Java/Others)
Total Submission(s): 6497 Accepted Submission(s): 1987

Problem Description
M斐波那契数列F[n]是一种整数数列，它的定义如下：

F[0] = a
F[1] = b
F[n] = F[n-1] * F[n-2] ( n > 1 )

现在给出a, b, n，你能求出F[n]的值吗？

Input
输入包含多组测试数据；
每组数据占一行，包含3个整数a, b, n（ 0 <= a, b, n <= 10^9 ）

Output
对每组测试数据请输出一个整数F[n]，由于F[n]可能很大，你只需输出F[n]对1000000007取模后的值即可，每组数据输出一行。

Sample Input
0 1 0
6 10 2

Sample Output
0
60

从题目我们可以很简单地就推出公式，a，b的指数是和菲波那切数列有关的，所以我们目标就变为快速求某一项的斐波那契值，再快速幂求答案。

对于这道题，a，b的指数非常大，所以我们当然不能直接对指数取模，为什么呢？
因为是错的，所以我们要用到欧拉降幂，来降低指数。

指数爆炸的时候就要降幂
就是求a^b mod c
可以转化为
a^(b mod phi©+phi©) mod c

phi就是欧拉函数，一般的求法为：

 ll phi(ll n)
{
     ll i,rea=n;
     for(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;
}

然而这道题我们可以不用求质数p的欧拉函数，为什么呢？
因为p是一个质数，我们可以从求欧拉函数的代码中可以看出，质数的欧拉函数就是本身减一。
从欧拉函数的定义我们也可以看出：对于一个数n，phi(n) 就是小于n的数当中与n互质的数的个数，因为p为质数，小于p的数肯定与p互质，所以每次取模时，我们对（p-1）取模就行。

AC代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll p=1e9+7;
ll a,b,n;
ll g[2][2],s[2][2],t[2][2];
void mul(ll a[][2],ll b[][2])
{
    memset(t,0,sizeof t);
    for(int k=0;k<2;k++)
        for(int i=0;i<2;i++)
            for(int j=0;j<2;j++)
                t[k][i]=((t[k][i]+a[k][j]*b[j][i]))%(p-1);
    for(int i=0;i<2;i++)    for(int j=0;j<2;j++)    a[i][j]=t[i][j];
}
ll qmi_e(ll n,ll a[][2])
{
    while(n)
    {
        if(n&1)    mul(a,g);
        n>>=1;
        mul(g,g);
    }
}
ll qmi(ll a,ll b)
{
    ll res=1;
    while(b)
    {
        if(b&1)    res=res*a%p;
        b>>=1;
        a=a*a%p;
    }
    return res;
}
int main()
{  
	while(cin>>a>>b>>n)
    {
        g[0][0]=1;g[0][1]=1;g[1][0]=1;g[1][1]=0;
        if(n==0)
        {
            cout<<a<<endl;continue;
        }
        if(n==1)
        {
            cout<<b<<endl;continue;
        }
        if(n==2)
        {
            cout<<a*b<<endl;continue;
        }
        for(int i=0;i<2;i++)    for(int j=0;j<2;j++)    s[i][j]=g[i][j];
        n-=3;
        qmi_e(n,s);
        ll k1=(s[0][0]+s[1][0])%(p-1);
        ll k2=s[0][0]%(p-1);
        ll res=qmi(a,k2)*qmi(b,k1)%p;
        cout<<res<<endl;
    }
    return 0;
}