zoukankan      html  css  js  c++  java
  • 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("
    ");
    		}printf("
    ");
    	}
    }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
    ",ans);
    	}return 0;
    }
    

      

  • 相关阅读:
    XNA入门教程(一)
    SQL透视表
    java 远程ftp建立文件夹
    费事数列——我的理解
    OOP
    OOP2
    河内之塔
    获取页面上TextBox并改变它的值
    RMAN学习之三:归档模式有备份,丢失控制文件。
    SQL Server Error: [DBNETLIB][ConnectionOpen (Connect()).]SQL Server 不存在或访问
  • 原文地址:https://www.cnblogs.com/joyouth/p/5381512.html
Copyright © 2011-2022 走看看