51nod 1195 斐波那契数列的循环节【斐波那契数列&&二次剩余&&欧拉判定准则】

解题思路：

先说明一下结论在下都不会证明，囧……。

对于一个正整数n，我们求Fib数模n的循环节的长度的方法如下：
$(1)将n分解质因数，即n=p_1^{k_1}p_2^{k_2}……p_m^{k_m}$
$(2)分别计算Fib数模p_i^{k_i}的循环节的长度，假设是x_1,x_2,……x_m$
$(3)则Fib数模n的循环节的长度为ans=lcm(x_1,x_2,……x_m)$
从上面三个步骤看来，貌似最困难的是第二步，那么我们如何求Fib数模$p^m$的循环节长度呢？

这里有一个优美的定理：

$Fib数模p^m的循环节长度等于G(p)*p^{m-1}，其中G(p)表示Fib数模p的循环节长度$

那么现在的关键在于如何求$G(p)$。

当$p\le5$，直接手算。

当$p>5$，又有结论如下：

$当5是模p的二次剩余，即方程x^2\equiv 5(mod$ $p)有解，那么循环节长度是(p-1)的约数，否则是2(p+1)的约数$

之所以与5有关是因为斐波那契通项公式里唯一的无理数是$\sqrt5，x$相当于起到了$\sqrt5$的作用。

接下来直接dfs枚举约数，加上矩阵快速幂判定即可，即满足$f[ans+1]=f[ans+2]=1$。注意矩阵乘法要手写卡常数。

---

这里提一下如何判定a是否是模质数p的二次剩余，即方程$x^2\equiv a(mod$ $p)$是否有解的判定方法。
该方法叫做欧拉判定准则：

$当a^{\frac{p-1}{2}}\equiv1(mod$ $p)，则a是模p的二次剩余，反之，即a^{\frac{p-1}{2}}\equiv-1(mod$ $p)则不是$

这个可以用费马小定理证明。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<cmath>
#include<ctime>
#define ll long long
using namespace std;

int getint()
{
    int i=0,f=1;char c;
    for(c=getchar();(c<'0'||c>'9')&&c!='-';c=getchar());
    if(c=='-')c=getchar(),f=-1;
    for(;c>='0'&&c<='9';c=getchar())i=(i<<3)+(i<<1)+c-'0';
    return i*f;
}

inline void W(long long x)
{
    x<0?putchar('-'),x=-x:0;
    static int buf[25];
    int tot=0;
    do{
        buf[tot++]=x%10,x/=10;
    }while(x);
    while(tot)putchar(buf[--tot]+'0');
}

const int N=100005;
const ll INF=1e18;
struct node
{
    int num,cnt;
    ll len;
}yz[50];
struct matrix
{
    int a[3][3];
    inline matrix I()
    {
        matrix res;
        res.a[1][1]=res.a[2][2]=1;
        res.a[1][2]=res.a[2][1]=0;
        return res;
    }
    inline friend matrix mul(const matrix &A,const matrix &B,const int &p)
    {
        matrix res;
        memset(res.a,0,sizeof(res.a));
        res.a[1][1]=(1ll*A.a[1][1]*B.a[1][1]+1ll*A.a[1][2]*B.a[2][1])%p;
        res.a[1][2]=(1ll*A.a[1][1]*B.a[1][2]+1ll*A.a[1][2]*B.a[2][2])%p;
        res.a[2][1]=(1ll*A.a[2][1]*B.a[1][1]+1ll*A.a[2][2]*B.a[2][1])%p;
        res.a[2][2]=(1ll*A.a[2][1]*B.a[1][2]+1ll*A.a[2][2]*B.a[2][2])%p;
        return res;
    }
    inline matrix Pow(int y,const int &p)
    {
        matrix res=I(),x=*this;
        for(;y;y>>=1,x=mul(x,x,p))
            if(y&1)res=mul(res,x,p);
        return res;
    }   
};
int T,n,tot;
ll t;
int pn,pri[N];
int len[N],tmp[N];
bool notp[N];

inline int ksm(int x,int y,int p)
{
    int res=1;
    for(;y;y>>=1,x=1ll*x*x%p)
        if(y&1)res=1ll*res*x%p;
    return res;
}

inline matrix Fib(int x,const int &p)
{
    matrix res,A;
    res.a[1][1]=res.a[2][1]=res.a[2][2]=0;
    res.a[1][2]=1;
    A.a[1][1]=0;
    A.a[1][2]=A.a[2][1]=A.a[2][2]=1;
    res=mul(res,A.Pow(x-1,p),p);
    return res;
}

inline void div(ll x)
{
    tot=0;
    for(int i=1;1ll*pri[i]*pri[i]<=x&&i<=pn;i++)
        if(x%pri[i]==0)
        {
            yz[++tot].num=pri[i],yz[tot].cnt=0,yz[tot].len=1;
            while(x%pri[i]==0)
            {
                x/=pri[i];
                yz[tot].cnt++;
            }
        }
    if(x!=1)yz[++tot].num=x,yz[tot].cnt=1,yz[tot].len=1;
}

void dfs(int x,ll num,int p)
{
    if(x==tot+1)
    {
        matrix A=Fib(num+2,p);
        if(A.a[1][1]==1&&A.a[1][2]==1)
            t=min(num,t);
        return;
    }
    dfs(x+1,num,p);
    ll d=1;
    for(int i=1;i<=yz[x].cnt;i++)
    {
        d*=yz[x].num;
        dfs(x+1,1ll*num*d,p);
    }
}

inline ll find(int p)
{
    if(p==2)return 3;
    if(p==3)return 8;
    if(p==5)return 20;
    t=INF; 
    if(ksm(5,(p-1)/2,p)==1)div(p-1);
    else div(2ll*(p+1));
    dfs(1,1,p);
    return t;
}

inline void sieve()
{
    int M=sqrt(1e9)+1;
    for(int i=2;i<=M;i++)
    {
        if(!notp[i])pri[++pn]=i,len[i]=find(i);
        for(int j=1;j<=pn;j++)
        {
            int k=i*pri[j];
            if(k>M)break;
            notp[k]=true;
            if(!i%pri[j])break;
        }
    }
}

ll gcd(ll a,ll b)
{
    return b?gcd(b,a%b):a;
}

inline void solve(int x)
{
    div(x);
    for(int i=1;i<=tot;i++)
    {
        for(int j=1;j<yz[i].cnt;j++)
            yz[i].len*=yz[i].num;
        if(1ll*yz[i].num*yz[i].num<=x)
            yz[i].len*=len[yz[i].num];
    }
    tmp[0]=tot;
    for(int i=1;i<=tot;i++)tmp[i]=yz[i].len;
    if(1ll*yz[tot].num*yz[tot].num>x)tmp[tmp[0]]=find(yz[tot].num);
    ll ans=1;
    for(int i=1;i<=tmp[0];i++)
        ans=1ll*ans*tmp[i]/gcd(ans,tmp[i]);
    W(ans),putchar('\n');
}

int main()
{
    //freopen("lx.in","r",stdin);
    //freopen("lx.out","w",stdout);
    sieve();
    T=getint();
    while(T--)
    {
        n=getint();
        solve(n);
    }
    return 0;
}