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

      

  • 相关阅读:
    每日一水 POJ8道水题
    编译和使用 MySQL C++ Connector
    j2ee model1模型完成分页逻辑的实现 详解!
    DB查询分析器访问EXCEL时,要在表名前后加上中括弧或双引号
    指向结构体变量的指针
    EOSS V3.0 企业运营支撑系统(基于RBAC原理的权限管理)
    MybatisGen1.0 Mybatis JavaBean Mapper生成工具
    The table name must be enclosed in double quotation marks or sqare bracket while accessing EXCEL by
    资源-Android:Android
    软件-开发软件:Android Studio
  • 原文地址:https://www.cnblogs.com/joyouth/p/5381512.html
Copyright © 2011-2022 走看看