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;
}