2016.4.11Test:problem3:龙卷风:Fibnacci的循环节,矩阵快速幂

原题HDU4794：Arnold

题目中(x,y)->((x+y)%n,(x+2*y)%n)=(x1,y1)，展开后是一个Fiboncci数列

所以矩阵转移回来等价于求Fiboncci数列每个模n意义下的循环节的LCM

对于一个正整数n，我们求Fib数模n的循环节的长度的方法如下：把n素因子分解，即n=p1^a1*p2^a2*...*pk^ak；分别计算Fib数模每个p^m的循环节长度，假设长度分别是x1,x2,...,xk；那么Fib模n的循环节长度L=lcm(x1,x2,...,xk)，Fib数模p^m的最小循环节长度等于G(p)*p^(m-1)，其中G(p)表示Fib数模素数p的最小循环节长度。可以看出我们现在最重要的就是求G(p)。我们利用如下定理：如果5是模p的二次剩余，那么循环节的的长度是p-1的因子，否则，循环节的长度是2(p+1)的因子。对于小于等于5的素数，我们直接特殊判断，loop(2)=3,loop(3)=8,loop(5)=20。我们可以先求出所有的因子，然后用矩阵快速幂来一个一个判断，这样时间复杂度不会很大。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#define LL unsigned long long
using namespace std;
const int maxb=1000010;
LL n,a[2][2],ret[2][2],sta[maxb],cnt=0,ret1,ret2;

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 powe(LL x,LL y,LL mod){
	LL ret=1;
	while (y){
		if (y&1) ret=ret*x%mod;
		x=x*x%mod; y>>=1;
	}return ret;
}

void Matrix_mul(LL a[2][2],LL b[2][2],LL c[2][2],LL mod){
	LL tmp[2][2]={0};
	for (int i=0;i<2;++i)
	    for (int j=0;j<2;++j)
			for (int k=0;k<2;++k)
			    tmp[i][j]=(tmp[i][j]+a[i][k]*b[k][j]%mod)%mod;
	for (int i=0;i<2;++i)
	    for (int j=0;j<2;++j)
	        c[i][j]=tmp[i][j];
}

void Matrix_powe(LL a[2][2],LL y,LL mod){
	memset(ret,0,sizeof(ret));
	ret[0][0]=ret[1][1]=1;
	while (y){
		if (y&1) Matrix_mul(ret,a,ret,mod);
		Matrix_mul(a,a,a,mod); y>>=1;
	}
}

void Fibonaci(LL p,LL mod){
	memset(a,0,sizeof(a));
	a[0][0]=1; a[0][1]=1;
	a[1][0]=1; a[1][1]=0;
	Matrix_powe(a,p-1,mod);
	ret1=(ret[1][0]+ret[1][1])%mod;
	ret2=(ret[0][0]+ret[0][1])%mod;
}

LL loop(LL mod){
	LL tmp; cnt=0;
	if (mod==2) return 3;
	if (mod==3) return 8;
	if (mod==5) return 20;
	if (powe(5,(mod-1)>>1,mod)==1) tmp=mod-1;
	else tmp=2*(mod+1);
	for (LL i=1;i*i<=tmp;++i)
	    if (tmp%i==0){
			LL x=i,y=tmp/i;
			Fibonaci(x,mod);
			if (ret1==0&&ret2==1) return x;
			if (y!=x) sta[cnt++]=y;
	    }
	while (cnt>0){
		Fibonaci(sta[--cnt],mod);
		if (ret1==0&&ret2==1) return sta[cnt];
	}return 0;
}

int main(){
	while (scanf("%lld",&n)!=EOF){
		LL tmp=n,ans=1;
		for (LL i=2;i*i<=tmp;++i)
		    if (tmp%i==0){
				LL len=loop(i),S=1;
				while(!(tmp%i))tmp/=i,S*=i;
				S/=i; S*=len;
				ans=lcm(ans,S);
		    }
		if (tmp>1){
			LL len=loop(tmp);
			ans=lcm(ans,len);
		}
		if (!(ans%2)) ans/=2;
		printf("%lld\n",ans);
	}
}