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

  1 #include<stdio.h>
  2 #include<string.h>
  3 #include<stdlib.h>
  4 #define N 1000005
  5 
  6 int visit[N],prime[N],K;
  7 long long P,Phi;
  8 
  9 struct matrix
 10 {
 11     long long A[2][2];
 12 };
 13 
 14 void intt()   // 找出1000000以内的素数
 15 {
 16     int i,j;
 17     memset(visit,0,sizeof(visit));
 18     for(i=2; i<=1000; i++)
 19     {
 20         if(visit[i]==0)
 21         {
 22             for(j=i+i; j<=1000000; j+=i)
 23             {
 24                 visit[j]=1;
 25             }
 26         }
 27     }
 28     K=0;
 29     for(j=2; j<=1000000; j++)
 30         if(visit[j]==0) prime[++K]=j;
 31 }
 32 
 33 matrix power(matrix ans1,matrix ans2)  // 矩阵乘法
 34 {
 35     int i,j,k;
 36     matrix ans;
 37     for(i=0; i<=1; i++)
 38     {
 39         for(j=0; j<=1; j++)
 40         {
 41             ans.A[i][j]=0;
 42             for(k=0; k<=1; k++)
 43             {
 44                 ans.A[i][j]+=ans1.A[i][k]*ans2.A[k][j];
 45                 if(ans.A[i][j]>Phi)
 46                 {
 47                     ans.A[i][j]=ans.A[i][j]%Phi+Phi;
 48                 }
 49             }
 50         }
 51     }
 52     return ans;
 53 }
 54 
 55 matrix mod(matrix un, long long k)  // 求矩阵的k次幂
 56 {
 57     matrix ans;
 58     ans.A[0][0]=1;
 59     ans.A[0][1]=0;
 60     ans.A[1][0]=0;
 61     ans.A[1][1]=1;
 62     while(k)
 63     {
 64         if(k%2) ans=power(ans,un);
 65         un=power(un,un);
 66         k/=2;
 67     }
 68     return ans;
 69 }
 70 
 71 long long mod1(long long a, long long k)  //求(a^k)%p
 72 {
 73     long long ans;
 74     ans=1;
 75     while(k)
 76     {
 77         if(k%2)
 78         {
 79             ans=ans*a;
 80             ans%=P;
 81         }
 82         a=a*a;
 83         a%=P;
 84         k/=2;
 85     }
 86     return ans%P;
 87 }
 88 
 89 int main()
 90 {
 91     int i,ncase,t;
 92     long long a,b,aa,bb,n,num,num1,num2;
 93     matrix init,ans;
 94 
 95     intt();
 96     scanf("%d",&ncase);
 97 
 98     for(t=1; t<=ncase; t++)
 99     {
100         scanf("%I64d%I64d%I64d%I64d",&a,&b,&P,&n);
101         printf("Case #%d: ",t);
102         if(n==1)
103         {
104             printf("%I64d\n",a%P);
105             continue;
106         }
107         else if(n==2)
108         {
109             printf("%I64d\n",b%P);
110             continue;
111         }
112         else if(n==3)
113         {
114             printf("%I64d\n",a*b%P);
115             continue;
116         }
117         if(P==1)
118         {
119             printf("0\n");
120             continue;
121         }
122 
123         // 初始化求斐波那契数的初始矩阵
124         init.A[0][0]=0;
125         init.A[0][1]=1;
126         init.A[1][0]=1;
127         init.A[1][1]=1;
128         //  A^B % C = A ^ ( B % phi[C] + phi[C] ) %C  ( B >= phi[C] ) ,phi[C]表示与C互质的数的个数
129         Phi=1;
130         num=P;
131 
132         for(i=1; i<=K; i++)
133         {
134             if(prime[i]>P) break;
135             if(P%prime[i]==0)
136             {
137                 Phi*=(prime[i]-1);
138                 num/=prime[i];
139             }
140         }
141         //phi[C]=C*(1-1/p1)*(1-1/p2)*...*(1-1/pt);p1,p2,...pt表示C的素因子
142         Phi*=num;//Phi表示phi[C] 在这里c = p
143 
144 
145         ans=mod(init,n-3);//求指数
146         num1=ans.A[1][1];//a的指数
147         num2=ans.A[0][1]+ans.A[1][1];//b的指数    求b的指数不是已经溢出了吗。？？？
148         if(num2>Phi) num2=num2%Phi+Phi;
149 
150         aa=mod1(a,num1);//a^num1%p;
151         bb=mod1(b,num2);//b^num2%p;
152 
153         printf("%I64d\n",aa*bb%P);
154     }
155     return 0;
156 }