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

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<cstdio>#include<cstring>#include<iostream>#include<algorithm>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;i<maxn;i++){        if(!vis[i]){            prime[s++]=i;        }        for(LL j=0;j<s&&i*prime[j]<maxn;j++){            vis[i*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;i<fac[step][1];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\n",solve(n));    }    return 0;}