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

解题思路：

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

对于一个正整数n，我们求Fib数模n的循环节的长度的方法如下：

(1)将n分解质因数，即n=pk11pk22……pkmm

(2)分别计算Fib数模pkii的循环节的长度，假设是x1,x2,……xm

(3)则Fib数模n的循环节的长度为ans=lcm(x1,x2,……xm)

从上面三个步骤看来，貌似最困难的是第二步，那么我们如何求Fib数模pm的循环节长度呢？

这里有一个优美的定理：

Fib数模pm的循环节长度等于G(p)∗pm−1，其中G(p)表示Fib数模p的循环节长度

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

当p≤5，直接手算。

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

当5是模p的二次剩余，即方程x2≡5(mod p)有解，那么循环节长度是(p−1)的约数，否则是2(p+1)的约数

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

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

这里提一下如何判定a是否是模质数p的二次剩余，即方程x2≡a(mod p)是否有解的判定方法。

该方法叫做欧拉判定准则：

当ap−12≡1(mod p)，则a是模p的二次剩余，反之，即ap−12≡−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
;
int len
,tmp
;
bool notp
;

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