HDU3978(再论斐波那契数列循环节)
题目：http://acm.hdu.edu.cn/showproblem.php?pid=3978
题意：给一个Fibnacci数列，所以有f(0)=f(1)=1,然后有一个函数：G(n,i)=f(G(n,i-1)),其中有G(n,0)=f(n),然后给
三个数n，k，p，求G(n,k)%p的值。
分析：对于本题其实思路还算简单，就是：我们把G(n,k)写成G(n,k)=f(f(f...f(n,0)))形式，然后一层一层找循环节。
当然这个问题很早以前就做过啊，详见这里：http://blog.csdn.net/acdreamers/article/details/8226873
所以我们关键是先一层一层找寻环节，然后倒着计算，我认为找寻环节似乎是最不好办的，不过上次对于这个问题我已经有详
细的说明：http://blog.csdn.net/acdreamers/article/details/10983813
再说一遍：如果对于一个数n，找fibnacci数列模n的循环节长度，就是先对n素因子分解，然后分别计算，最后按照最小公倍
数合并，所以主要是求fibnacci模素数p的循环节长度了。这个，上次我说了有一个优美的定理。
定理：对于大于5的素数p，如果5是模p的二次剩余，那么fibnacci模p的循环节长度是p-1的因子，否则，如果5是模p的二次
非剩余，那么fibnacci数列模p的循环节长度是2*（p+1）的因子。
关于这个定理的证明就直接用费马小定理，fibnacci数列的公式，二次剩余和二次域证明即可。
那么，对于小于5的素数我们特殊处理即可，loop(2)=3,loop(3)=8,loop(5)=20
#include <iostream>#include <string.h>#include <algorithm>#include <stdio.h>#include <math.h>using namespace std;typedef long long LL;const int M=2;struct Matrix{ LL m[M][M];};Matrix per= {1,0,0,1};Matrix multi(Matrix a,Matrix b,LL MOD){ Matrix c; int i,j,k; for(i=0; i<M; i++) { for(j=0; j<M; j++) { c.m[i][j]=0; for(k=0; k<M; k++) { c.m[i][j]+=a.m[i][k]*b.m[k][j]%MOD; } c.m[i][j]%=MOD; } } return c;}Matrix power(Matrix a,LL k,LL MOD){ Matrix ans=per,p=a; while(k) { if(k&1) { ans=multi(ans,p,MOD); k--; } k>>=1; p=multi(p,p,MOD); } return ans;}LL gcd(LL a,LL b){ return b? gcd(b,a%b):a;}LL quick_mod(LL a,LL b,LL m){ LL ans=1; a%=m; while(b) { if(b&1) { ans=ans*a%m; b--; } b>>=1; a=a*a%m; } return ans;}//勒让德符号int legendre(int a,int p){ if(quick_mod(a,(p-1)>>1,p)==1) return 1; else return -1;}const int N=1000005;const int NN=50005;bool prime[N];int p[N];int num[NN],pri[NN];int num1[NN],pri1[NN];int arr[NN];int loop[N];int k,cnt,c;void isprime(){ k=0; int i,j; memset(prime,true,sizeof(prime)); for(i=2; i<N; i++) { if(prime[i]) { p[k++]=i; for(j=i+i; j<N; j+=i) { prime[j]=false; } } }}void find(int n,int pri[],int num[]){ cnt=0; int t=(int)sqrt(1.0*n); for(int i=0; p[i]<=t; i++) { if(n%p[i]==0) { int a=0; pri[cnt]=p[i]; while(n%p[i]==0) { a++; n/=p[i]; } num[cnt]=a; cnt++; } } if(n>1) { pri[cnt]=n; num[cnt]=1; cnt++; }}void dfs(int dept,int product=1){ if(dept==cnt) { arr[c++]=product; return; } for(int i=0; i<=num1[dept]; i++) { dfs(dept+1,product); product*=pri1[dept]; }}int find_loop(int n){ find(n,pri,num); int cnt1=cnt; LL ans=1; for(int i=0; i<cnt1; i++) { c=0; int record=1; if(pri[i]==2) record=3; else if(pri[i]==3) record=8; else if(pri[i]==5) record=20; else { if(legendre(5,pri[i])==1) find(pri[i]-1,pri1,num1); else find(2*(pri[i]+1),pri1,num1); dfs(0,1); sort(arr,arr+c); for(int k=0; k<c; k++) { Matrix A; A.m[0][0]=1; A.m[0][1]=1; A.m[1][0]=1; A.m[1][1]=0; Matrix a=power(A,arr[k]-1,pri[i]); int x=(a.m[0][0]+a.m[0][1])%pri[i]; int y=(a.m[1][0]+a.m[1][1])%pri[i]; if(x==1&&y==0) { record=arr[k]; break; } } } for(int k=1; k<num[i]; k++) record*=pri[i]; ans=ans/gcd(ans,record)*record; } return ans;}void Solve(int p,int k){ loop[0]=p; for(int i=1; i<=k; i++) loop[i]=find_loop(loop[i-1]);}int work(int n,int k,int p){ int t=n; LL ret,MOD; Matrix ans; Matrix A; A.m[0][0]=1; A.m[0][1]=1; A.m[1][0]=1; A.m[1][1]=0; Solve(p,k); for(int i=k; i>=0; i--) { MOD=loop[i]; ans=power(A,t,MOD); ret=(ans.m[1][0]+ans.m[1][1])%MOD; t=ret; } return ret;}int main(){ isprime(); int T,n,k,p,tt=1; scanf("%d",&T); while(T--) { scanf("%d%d%d",&n,&k,&p); printf("Case #%d: %d\n",tt++,work(n,k,p)); } return 0;}