HDU 3221 矩阵快速幂+欧拉函数+降幂公式降幂

装载自：http://www.cnblogs.com/183zyz/archive/2012/05/11/2495401.html 题目让求一个函数调用了多少次。公式比较好推。f[n] = f[n-1]*f[n-2]。然后a和b系数都是呈斐波那契规律增长的。需要先保存下来指数。但是太大了。在这里不能用小费马定理。要用降幂公式取模、
(A^x)%C=A^(x%phi(C)+phi(C))%C（x>=phi(C)） Phi[C]表示不大于C的数中与C互质的数的个数，可以用欧拉函数来求。

矩阵快速幂也不熟、。觉得很难。

#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#define N 1000005

int visit[N],prime[N],K;
long long P,Phi;

struct matrix
{
    long long A[2][2];
};

void intt()   // 找出1000000以内的素数
{
    int i,j;
    memset(visit,0,sizeof(visit));
    for(i=2; i<=1000; i++)
    {
        if(visit[i]==0)
        {
            for(j=i+i; j<=1000000; j+=i)
            {
                visit[j]=1;
            }
        }
    }
    K=0;
    for(j=2; j<=1000000; j++)
        if(visit[j]==0) prime[++K]=j;
}

matrix power(matrix ans1,matrix ans2)  // 矩阵乘法
{
    int i,j,k;
    matrix ans;
    for(i=0; i<=1; i++)
    {
        for(j=0; j<=1; j++)
        {
            ans.A[i][j]=0;
            for(k=0; k<=1; k++)
            {
                ans.A[i][j]+=ans1.A[i][k]*ans2.A[k][j];
                if(ans.A[i][j]>Phi)
                {
                    ans.A[i][j]=ans.A[i][j]%Phi+Phi;
                }
            }
        }
    }
    return ans;
}

matrix mod(matrix un, long long k)  // 求矩阵的k次幂
{
    matrix ans;
    ans.A[0][0]=1;
    ans.A[0][1]=0;
    ans.A[1][0]=0;
    ans.A[1][1]=1;
    while(k)
    {
        if(k%2) ans=power(ans,un);
        un=power(un,un);
        k/=2;
    }
    return ans;
}

long long mod1(long long a, long long k)  //求(a^k)%p
{
    long long ans;
    ans=1;
    while(k)
    {
        if(k%2)
        {
            ans=ans*a;
            ans%=P;
        }
        a=a*a;
        a%=P;
        k/=2;
    }
    return ans%P;
}

int main()
{
    int i,ncase,t;
    long long a,b,aa,bb,n,num,num1,num2;
    matrix init,ans;

    intt();
    scanf("%d",&ncase);

    for(t=1; t<=ncase; t++)
    {
        scanf("%I64d%I64d%I64d%I64d",&a,&b,&P,&n);
        printf("Case #%d: ",t);
        if(n==1)
        {
            printf("%I64d
",a%P);
            continue;
        }
        else if(n==2)
        {
            printf("%I64d
",b%P);
            continue;
        }
        else if(n==3)
        {
            printf("%I64d
",a*b%P);
            continue;
        }
        if(P==1)
        {
            printf("0
");
            continue;
        }

        // 初始化求斐波那契数的初始矩阵
        init.A[0][0]=0;
        init.A[0][1]=1;
        init.A[1][0]=1;
        init.A[1][1]=1;
        //  A^B % C = A ^ ( B % phi[C] + phi[C] ) %C  ( B >= phi[C] ) ,phi[C]表示与C互质的数的个数
        Phi=1;
        num=P;

        for(i=1; i<=K; i++)
        {
            if(prime[i]>P) break;
            if(P%prime[i]==0)
            {
                Phi*=(prime[i]-1);
                num/=prime[i];
            }
        }
        //phi[C]=C*(1-1/p1)*(1-1/p2)*...*(1-1/pt);p1,p2,...pt表示C的素因子
        Phi*=num;//Phi表示phi[C] 在这里c = p


        ans=mod(init,n-3);//求指数
        num1=ans.A[1][1];//a的指数
        num2=ans.A[0][1]+ans.A[1][1];//b的指数    求b的指数不是已经溢出了吗。？？？
        if(num2>Phi) num2=num2%Phi+Phi;

        aa=mod1(a,num1);//a^num1%p;
        bb=mod1(b,num2);//b^num2%p;

        printf("%I64d
",aa*bb%P);
    }
    return 0;
}