hdu 4794 FIb求循环节

很容易看出来这道题是求模n意义下fib数列的最小循环节

对于fib数列的最小循环节的求法，我们可以这样：

1、令n=p1^m1 * p2^m2 * p3^m3……

2、分别计算fib数列在模p1^m1,p2^m2……意义下的最小循环节

3、模n意义下的最小循环节为2步骤各循环节的LCM

首先步骤三是非常容易证明的，CRT直接可以看出来

我们考虑如何计算p1^m1的循环节

这里有一个定理:G(p1^m1)=G(p1)*p1^(m-1)

我们只需要计算G(p)就可以了

我们知道对于Fib数列，满足f(i)=1.0/sqrt(5) *( (1+sqrt(5))/2)^i - ((1-sqrt(5))/2)^i）

那么在模p意义下存在以下两种情况：

1、5是模p意义下的二次剩余，即满足欧拉判别 5^((p-1)/2)=1

那么原式中的sqrt(5)在模意义下可以用整数x替代

由费马小定理很容易知道F(p-1)=0

又知道

((1+x)/2)^p=(1+x)/2

((1-x)/2)^p=(1-x)/2

两式做差等于x

所以F(p)=1

因此我们知道p-1是Fib数列的一个循环节，那么其最小循环节一定是p-1的约数

枚举+矩阵乘法判断即可

2、5不是模p意义下的二次剩余，则存在5^((p-1)/2)=-1

那么由二项式定理得

(1+sqrt(5))^p=1+(sqrt(5))^p

(1-sqrt(5))^p=1-(sqrt(5))^p

(因为其系数的组合数里含有p这个因子的时候值为0）

两式做差得2*sqrt(5)^p

进而得F(p)=(sqrt(5)/2)^(p-1)

我们知道在模p意义下2^(p-1)的逆元为1

所以F(p)=(sqrt(5))^(p-1)=5^((p-1)/2)=-1

同理，我们可以得到F(p+1)=0

设Fib数列在模意义下第一个F(i)=0且i不等于0的位置为i

设F(i-1)=a 易证(a,p)=1

则一定有i是p+1的约数

且F(p)=a^j*F(p-i*j)

显然F(p-i*j)=a，所以F(p)=a^(j+1）

进而知道F(2*p+1)=a^(2*j+2)=(-1)^2=1

F(2*p+2)=0

所以2*(p+1)是Fib数列的一个循环节，那么其最小循环节一定是2*(p+1)的约数

枚举+矩阵乘法判断即可

 

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<iostream>
#include<algorithm>
using namespace std;

const int maxn=20010;
typedef long long LL;
LL n,mod,now;
LL a,b,ans;
LL st[maxn],top=0;
struct Matrix{
	LL a[2][2];
	void clear(){memset(a,0,sizeof(a));}
	void print(){
		for(int i=0;i<2;++i){
			for(int j=0;j<2;++j){
				printf("%lld ",a[i][j]);
			}printf("\n");
		}printf("\n");
	}
}A;
LL mul(LL a,LL b){
	LL s=0;
	while(b){
		if(b&1)s=(s+a)%mod;
		a=(a<<1)%mod;b>>=1;
	}return s;
}
LL GCD(LL a,LL b){return b==0?a:GCD(b,a%b);}
LL LCM(LL a,LL b){return a*b/GCD(a,b);}
LL Get_pow(LL v,LL p){
	LL tmp=1;
	while(p){
		if(p&1)tmp=mul(tmp,v);
		v=mul(v,v);p>>=1;
	}return tmp;
}
Matrix operator *(const Matrix &A,const Matrix &B){
	Matrix C;C.clear();
	for(int i=0;i<2;++i){
		for(int j=0;j<2;++j){
			for(int k=0;k<2;++k){
				C.a[i][j]=C.a[i][j]+mul(A.a[i][k],B.a[k][j]);
				if(C.a[i][j]>=mod)C.a[i][j]-=mod;
			}
		}
	}return C;
}
Matrix pow_mod(Matrix v,LL p){
	Matrix tmp;tmp.clear();
	tmp.a[0][0]=tmp.a[1][1]=1;
	while(p){
		if(p&1)tmp=tmp*v;
		v=v*v;p>>=1;
	}return tmp;
}
void F(LL p){
	A.a[0][0]=1;A.a[0][1]=1;
	A.a[1][0]=1;A.a[1][1]=0;
	A=pow_mod(A,p-1);
	b=(A.a[0][0]+A.a[0][1])%mod;
	a=(A.a[1][0]+A.a[1][1])%mod;
}
LL Get_loop(LL k){
	if(k==2)return 3;
	else if(k==3)return 8;
	else if(k==5)return 20;
	LL p;mod=k;top=0;
	if(Get_pow(5LL,(mod-1)>>1)==1)p=mod-1;
	else p=2*(mod+1);
	for(LL i=1;i*i<=p;++i){
		if(p%i==0){
			LL x=i,y=p/i;
			F(x);
			if(a==0&&b==1)return x;
			if(y!=x)st[++top]=y;
		}
	}
	while(top){
		F(st[top]);
		if(a==0&&b==1)return st[top];
		top--;
	}return 0;
}

int main(){
	while(scanf("%lld",&n)==1){
		now=n;ans=1;
		for(LL i=2;i*i<=now;++i){
			if(now%i==0){
				LL S=1;
				while(now%i==0)now/=i,S*=i;
				S/=i;S=S*Get_loop(i);
				ans=LCM(ans,S);
			}
		}
		if(now>1){
			LL S=Get_loop(now);
			ans=LCM(ans,S);
		}
		if(ans%2==0)ans/=2;
		printf("%lld\n",ans);
	}return 0;
}