HDU 3221Brute-force Algorithm(降幂公式 神似hdu4549)

Brute-force AlgorithmTime Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 1940    Accepted Submission(s): 475


Problem DescriptionProfessor 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. 
InputThere are multiple test cases. The first line of the input contains an integer T, meaning the number of the test cases.

For each test cases, there are four integers a, b, P and n in a single line.
You can assume that 1≤n≤1000000000, 1≤P≤1000000, 0≤a, b<1000000. 
OutputFor each test case, output the answer with case number in a single line. 
Sample Input

33 4 10 34 5 13 53 2 19 100

 
Sample Output

Case #1: 2Case #2: 11Case #3: 12

 
Source2009 Asia Shanghai Regional Contest Host by DHU 

题目大意：这个题目主要开始感觉意思难懂一点，实际上就是迭代的思想，f(1)=a,  f(2)=b,  f(3)=f(2,f(1))=f(1)*f(2)=a*b,  f(4)=f(3,f(2))=f(3)*f(2)=a*b^2,  f(5)=a^2*b^3.....后面的规律就一目了然了吧，然后给你一个地址。hdu4549M斐波那契数列题目真的是惊人的相似.
      解题思路：只是两者有一点不同，前面那个hdu4549直接%1e9+7，而小费马定理a^(p-1)%p==1,(a,p互质，且p为质数)，很显然p==1e9+7满足这个条件，所以可以用小费马定理降幂：a^t%p==a^(t%(p-1))。可以理解为p-1为一个循环节。 但是这个题目就不是这样了，这也是两点的不同之处！这个p可能不是质数，因此费马小定理失效，但是可以用降幂公式：(A^x)%C=A^(x%phi(C)+phi(C))%C（x>=phi(C)）.主要思想就是用矩阵的快速幂把a,b的指数算出来，注意算指数的时候结合降幂公式，然后最后直接用快速幂算出结果即可。
      题目地址：Brute-force Algorithm
如果对矩阵的快速幂那一块或者开始推出fibo数列这块不清楚的，可以看HDU 4549M斐波那契数列(矩阵快速幂+费马小定理)

AC代码：

#include<iostream>#include<cstring>#include<string>#include<cstdio>#include<cmath>using namespace std;__int64 mo,phi;__int64 ret[2][2],tmp[2][2],p[2][2];__int64 n;void init()  //初始化{    ret[0][0]=1; ret[0][1]=1;    ret[1][0]=1; ret[1][1]=0;    p[0][0]=1; p[0][1]=1;    p[1][0]=1; p[1][1]=0;}void cal1()  //!(n&1){    int i,j,k;    for(i=0;i<2;i++)        for(j=0;j<2;j++)        {            tmp[i][j]=p[i][j];            p[i][j]=0;        }    for(i=0;i<2;i++)        for(j=0;j<2;j++)          for(k=0;k<2;k++)          {                p[i][j]=p[i][j]+tmp[i][k]*tmp[k][j];                if(p[i][j]>=phi)                    p[i][j]=p[i][j]%phi+phi;          }}void cal2()  //n&1{    int i,j,k;    for(i=0;i<2;i++)        for(j=0;j<2;j++)        {            tmp[i][j]=ret[i][j];            ret[i][j]=0;        }    for(i=0;i<2;i++)        for(j=0;j<2;j++)          for(k=0;k<2;k++)          {              ret[i][j]=ret[i][j]+tmp[i][k]*p[k][j];              if(ret[i][j]>=phi)                ret[i][j]=ret[i][j]%phi+phi;          }}void fastmi()  //矩阵的快速幂{    init();    n-=3;    while(n)    {        if(n&1)            cal2();        cal1();        n>>=1;    }}__int64 pow(__int64 base,__int64 p)  //快速幂{    __int64 ans=1;    while(p)    {        if(p&1)            ans=(ans*base)%mo;        base=(base*base)%mo;        p>>=1;    }    return ans;}__int64 geteuler(__int64 n){    __int64 m=sqrt(n+0.5),ans=n,i;    for(i=2;i<=m;i++)        if(n%i==0)        {            ans=ans/i*(i-1);            while(n%i==0)                n/=i;        }    if(n>1)        ans=ans/n*(n-1);    return ans;}int main(){    __int64 a,b;    int i,tes;    scanf("%d",&tes);    for(i=1;i<=tes;i++)    {        scanf("%I64d%I64d%I64d%I64d",&a,&b,&mo,&n);        n--;        phi=geteuler(mo);   //得到欧拉值        __int64 ans1,ans2,res1,res2;        if(n==0) ans1=1,ans2=0;        else if(n==1) ans1=0,ans2=1;        else if(n==2) ans1=1,ans2=1;        else        {            fastmi();            ans2=(ret[0][0]+ret[0][1]);  //b的次数            ans1=(ret[1][0]+ret[1][1]);  //a的次数        }        //printf("%I64d %I64d\n",ans1,ans2);        if(ans1>=phi)            ans1=ans1%phi+phi;        if(ans2>=phi)            ans2=ans2%phi+phi;        res1=pow(a,ans1);        res2=pow(b,ans2);        __int64 res=(res1*res2)%mo;        printf("Case #%d: %I64d\n",i,res);    }    return 0;}//0MS