51 nod 1195 斐波那契数列的循环节

题目来源： Spoj 基准时间限制：1 秒 空间限制：131072 KB 分值: 640 难度：8级算法题 斐波那契数列Mod 一个数N会得到一个新的数列，根据同余可以得知，这个数列中的数会出现循环。例如：
F (mod 4) = 0 1 1 2 3 1 …
F (mod 5) = 0 1 1 2 3 0 3 3 1 4 0 4 4 3 2 0 2 2 4 1 … 后面的数会出现循环，因此Fib Mod 4的循环节长度是6，Fib Mod 5的循环节长度是20。
给出一个数N，求斐波那契数列Mod N的循环节的长度。 Input 第1行：一个数T，表示后面用作输入测试的数的数量。（1 <= T <= 10000) 第2 - T + 1行：每行1个数N。（1 <= N <= 10^9) Output 共T行，每行一个数，对应循环节的长度。 Input示例 2
4
5 Output示例 6
20 按照唐老师题解写的，交了20几次，一直出现超时的问题，这题写了几天，断断续续的了解知识点和解决复杂度的问题。
以下是唐老师题解：/*
建议先了解一下二次剩余方面的知识，再读一下这篇论文，其中讲了一些关于斐波那契数列模意义下循环节的性质。
Charles W. Campbell II. The Period of the Fibonacci Sequence Modulo j. 2007. 主要用到几个性质：
1. 对于与5互质的质数p，如果5是模p的二次剩余，那么模p意义下的循环节长度是(p-1)的因子。
2. 对于与5互质的质数p，如果5是模p的二次非剩余，那么模p意义下的循环节长度是(2p+2)的因子。
3. 对于模质数的幂p^k意义下的循环节，其值为模p意义下的循环节长度乘p^(k-1)。
4. 对于模x意义下的循环节，如果x被质因数分解为p1^k1*p2^k2*…*pm^km，则循环节是模每个质数的幂意义下的循环节的最小公倍数。 所以算法流程就是将N分解为质数的幂的乘积，分别求解，对于模5的循环节可以预先算出，因为它必然是20的因子。
对于每个模质数p的情况，枚举循环节的倍数（(p-1)或(2p+2)）的因子，计算出最小循环节。 需要一些技巧将时间复杂度简化。
1. 将一些取模的操作用加减法代替；输入输出优化。
2. 筛出一些质数，加快质因数分解的操作，例如分解n的复杂度可以由O(sqrt(n))变成O(sqrt(n)/ln(n))。
3. 算模质数p意义下的循环节时，由于斐波那契数列的”循环节的倍数”次项在模意义下相等，所以可以用类似计算欧拉函数的形式将时间复杂度由O(sqrt(n))变成O(log(n))。
4. 斐波那契数列是线性递推数列，利用2*2的矩阵进行快速幂，每次的矩阵乘法是8次加法、8次乘法，但是使用二元组线性表示后，每次的二元组乘法是3次加法、4次乘法，可以节省一半时间，采用的公式可以看作是Fib(n+m)=Fib(n-1)*Fib(m)+Fib(n)*Fib(m+1)。
5. 利用记忆化搜索技巧减少小质数情况的计算。
其中，使用2、3项就已经可以在允许的时间范围内得到正确答案了，不过结合1、2、4、5也可以卡常数蹭过（主要是4）。*/
代码的注释 我就不多加了，唐老师的算法很清楚了，代码自己理解吧。 #include #include #include #include using namespace std; typedef long long LL; const LL maxn=1e5+5; LL dp[maxn*10]; LL prime[maxn],s=0; bool vis[maxn]; void init_prime(){ for(LL i=2;iif(!vis[i]){ prime[s++]=i; } for(LL j=0;j<s&&i*prime[j]*prime[j]]=1; if(i%prime[j]==0){ break; } } } return ; } LL pow_mod(LL a1,LL b1){ LL ans=1; while(b1){ if(b1&1){ ans=ans*a1; } b1>>=1; a1*=a1; } return ans; } LL pow_mod2(LL a,LL b,LL p){ LL ans=1; while(b){ if(b&1){ ans=ans*a%p; } b>>=1; a=a*a%p; } return ans; } LL gcd(LL a,LL b){ return b?gcd(b,a%b):a; } bool f(LL n,LL p){ if(pow_mod2(n,(p-1)>>1,p)!=1){ return false; } return true; } struct matrix{ LL x1,x2,x3,x4; }; matrix matrix_a,matrix_b; matrix M2(matrix aa,matrix bb,LL mod){ matrix tmp; tmp.x1=(aa.x1*bb.x1%mod+aa.x2*bb.x3%mod)%mod; tmp.x2=(aa.x1*bb.x2%mod+aa.x2*bb.x4%mod)%mod; tmp.x3=(aa.x3*bb.x1%mod+aa.x4*bb.x3%mod)%mod; tmp.x4=(aa.x3*bb.x2%mod+aa.x4*bb.x4%mod)%mod; return tmp; } matrix M(LL n,LL mod){ matrix a,b; a=matrix_a;b=matrix_b; while(n){ if(n&1){ b=M2(b,a,mod); } n>>=1; a=M2(a,a,mod); } return b; } /*bool g(matrix t){ return t.m[0][0]==1&&t.m[0][1]==0&&t.m[1][0]==0&&t.m[1][1]==1; }*/ LL fac[100][2],l,x,fs[1000]; void dfs(LL count,LL step){ if(step==l){ fs[x++]=count; return ; } LL sum=1; for(LL i=0;i1];i++){ sum*=fac[step][0]; dfs(count*sum,step+1); } dfs(count,step+1); } LL solve2(LL p){ if(p<1e6&&dp[p]){ return dp[p]; } bool ok=f(5,p);//判断5是否为p的二次剩余 LL t; if(ok){ t=p-1; } else{ t=2*p+2; } l=0; for(LL i=0;i<s;i++){ if(prime[i]>t/prime[i]){ break; } if(t%prime[i]==0){ LL count=0; fac[l][0]=prime[i]; while(t%prime[i]==0){ count++;t/=prime[i]; } fac[l++][1]=count; } } if(t>1){ fac[l][0]=t; fac[l++][1]=1; } x=0; dfs(1,0);//求t的因子 sort(fs,fs+x); for(LL i=0;i<x;i++){ matrix m1=M(fs[i],p); if(m1.x1==m1.x4&&m1.x1==1&&m1.x2==m1.x3&&m1.x2==0){ if(p<1e6){ dp[p]=fs[i]; } return fs[i]; } } } LL solve(LL n){ LL ans=1,cnt; for(LL i=0;i<s;i++){ if(prime[i]>n/prime[i]){ break; } if(n%prime[i]==0){ LL count=0; while(n%prime[i]==0){ count++;n/=prime[i]; } cnt=pow_mod(prime[i],count-1); cnt*=solve2(prime[i]); ans=(ans/gcd(ans,cnt))*cnt; } } if(n>1){ cnt=1; cnt*=solve2(n); ans=ans/gcd(ans,cnt)*cnt; } return ans; } int main(){ LL t,n; init_prime(); matrix_a.x1=matrix_a.x2=matrix_a.x3=1; matrix_a.x4=0; matrix_b.x1=matrix_b.x4=1; matrix_b.x2=matrix_b.x3=0; dp[2]=3;dp[3]=8;dp[5]=20; scanf("%I64d",&t); while(t--){ scanf("%I64d",&n); printf("%I64d ",solve(n)); } return 0; }