HDU 3977 Evil teacher（斐波那契数列循环节）

Description
求斐波那契数列模p的最小循环节m(m>0)
Input
第一行为一整数T表示用例组数，每组用例为一整数p(T<=20,1<=p<=2*10^9,p的最大素因子不超过10^6)
Output
求斐波那契数列模p的最小循环节m
Sample Input
5
11
19
61
17
67890
Sample Output
Case #1: 10
Case #2: 18
Case #3: 60
Case #4: 36
Case #5: 4440
Solution
首先我们知道斐波那契数列作为广义斐波那契数列f(n)=a*f(n-1)+b*f(n-2)的一个特殊情况在模素数p的最小循环节len满足以下结论（其中c=a^2+4*b^2=5）
1.当c是模p的二次剩余时，最小循环节是p-1的因子
2.当c不是模p的二次剩余时，最小循环节是(p-1)(p+1)的因子
3.当c=p时，结论不定，此处暴力找一下可以知道p=5时最小循环节是20
但是此处的p不是素数，先对p素因子分解，p=p1^k1*p2^k2*….*pm^km，对于每个素因子pi可以采用上面的结论求出在模pi下的最小循环节len[pi]，而模pi^ki的最小循环节在pi<=10^6的情况下经检验是len[pi]*pi^(ki-1)，那么对于每个pi^ki就可以求出其最小循环节，求出的若干最小循环节的lcm即为答案（可以用中国剩余定理验证）
Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int factor[33][2],res;
ll len[1111111];
ll gcd(ll a,ll b)
{
    return b?gcd(b,a%b):a; 
}
ll lcm(ll a,ll b)
{
    return a/gcd(a,b)*b;
}
struct Mat
{
    ll mat[3][3];//矩阵 
    int row,col;//矩阵行列数 
};
Mat mod_mul(Mat a,Mat b,ll p)//矩阵乘法 
{
    Mat ans;
    ans.row=a.row;
    ans.col=b.col;
    memset(ans.mat,0,sizeof(ans.mat));
    for(int i=0;i<ans.row;i++)      
        for(int k=0;k<a.col;k++)
            if(a.mat[i][k])
                for(int j=0;j<ans.col;j++)
                {
                    ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                    ans.mat[i][j]%=p;
                }
    return ans;
}
Mat Mat_mod_pow(Mat a,ll k,ll p)//矩阵快速幂 
{
    Mat ans;
    ans.row=a.row;
    ans.col=a.col;
    for(int i=0;i<a.row;i++)
        for(int j=0;j<a.col;j++)
            ans.mat[i][j]=(i==j);
    while(k)
    {
        if(k&1)ans=mod_mul(ans,a,p);
        a=mod_mul(a,a,p);
        k>>=1;
    }
    return ans;
}
ll deal(int x)
{
    if(len[x]!=-1)return len[x];
    Mat A;
    A.mat[0][0]=A.mat[0][1]=A.mat[1][0]=1,A.mat[1][1]=0;
    A.row=A.col=2;
    ll p=1ll*(x-1)*(x+1);
    ll ans=p;
    for(ll i=2;i*i<=p;i++)
        if(p%i==0)
        {
            Mat B=Mat_mod_pow(A,i,x);
            if(B.mat[0][0]==1&&B.mat[0][1]==0&&B.mat[1][0]==0&&B.mat[1][1]==1)
                ans=min(ans,i);
            B=Mat_mod_pow(A,p/i,x);
            if(B.mat[0][0]==1&&B.mat[0][1]==0&&B.mat[1][0]==0&&B.mat[1][1]==1)
                ans=min(ans,p/i);
        }
    return len[x]=ans;
}
void Get_Factor(int p)
{
    res=0;
    for(int i=2;i*i<=p;i++)
        if(p%i==0)
        {
            factor[res][0]=i,factor[res][1]=0;
            while(p%i==0)p/=i,factor[res][1]++;
            res++;
        }
    if(p>1)factor[res][0]=p,factor[res++][1]=1;
}
int main()
{
    memset(len,-1,sizeof(len));
    len[1]=1,len[2]=3,len[3]=8,len[4]=6,len[5]=20;
    int T,p,Case=1;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d",&p);
        Get_Factor(p);
        ll ans=1;
        for(int i=0;i<res;i++)
        {
            ll temp=deal(factor[i][0]);
            for(int j=1;j<factor[i][1];j++)temp*=factor[i][0];
            ans=lcm(ans,temp);
        }
        printf("Case #%d: %I64d\n",Case++,ans);
    }
    return 0;
}