hdu 3221 Brute-force Algorithm

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3221

矩阵乘法和数学的结合。

由题目意思可知：f
=f[n-1]*f[n-2];

f的前面几项可以罗列出来：

a^1*b^0,a^0*b^1,a^1*b^1,a^1*b^2,a^2*b^3....

可以发现a的指数和b的指数均类似于斐波那契数列。

用矩阵的快速幂可以很快的求出第n项a和b的指数分别是多少。

但是这个指数会非常大，存不下来，需要对一个数去模。

这里需要用到一个公式:

A^B%C=A^( B%Phi[C] + Phi[C] )%C (B>=Phi[C])

Phi[C]表示不大于C的数中与C互质的数的个数，可以用欧拉函数来求：

找到C的所有素因子。

Phi[C]=C*(1-1/q1)*(1-1/q2)*(1-1/q3)*....*(1-1-qk);

q1,q2,q3...qk表示C的素因子。

到这里基本上就能解决了。

code:

View Code

# include<stdio.h>
# include<string.h>
# include<stdlib.h>
# define N 1000005
int visit
,prime
,K;
__int64 P,Phi;
struct matrix{
__int64 A[2][2];
};
void intt()
{
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,__int64 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;
}
__int64 mod1(__int64 a,__int64 k)
{
__int64 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;
__int64 a,b,aa,bb,n,num,num1,num2;
//freopen("3221.in","r",stdin);
//freopen("32211.out","w",stdout);
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\n",a%P);continue;}
else if(n==2) {printf("%I64d\n",b%P); continue;}
else if(n==3) {printf("%I64d\n",a*b%P);continue;}
if(P==1) {printf("0\n");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]

ans=mod(init,n-3);///求指数

num1=ans.A[1][1];///a的指数
num2=ans.A[0][1]+ans.A[1][1];///b的指数
if(num2>Phi) num2=num2%Phi+Phi;

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

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