[斐波那契循环节 数学技巧] HDU 3977 Evil teacher

找到了两种做法 不同的关键在于对于一个质数p的循环节的不同求法

法一：

来自http://blog.csdn.net/prime7/article/details/11017111

分析过程：首先我们知道fib数列模p如果出现了连续的1，0就意味这着开始循环了，因为接下来的项就是1 1 2 3 5等等。

那么很显然如果在第k位第一次出现了1，0，那么对于以后的1，0都可以表示为k*m。

 

那么，现在我们考虑如果fib数列模p在第pos位第一次出现了0，那么设0前面的那个数为a，则接下来的序列将是a，0，a，

a，2a，3a，5a，8a，....。可以看出a的系数就是一个fib数列，那么我们就可以得到fib(k+i)%p=a*fib(i)%p，其中i满

足0<i<k，所以进一步可以得到fib(i)=[a^j*fib(i-k*j)]%p。

那么我们现在先来说说如何求fib数模一个正整数n的循环节长度：

对于这个问题，我们先对n进行素因子分解，得到：，然后先对每一个形如p^k的数计算循环节，然后它们

的最小公倍数就是n的循环节长度（当然这里面涉及到CRT等等方面的基础）。那么现在问题就是计算p^k的循环节，这个问题

可以进一步简化成求G(p)*p^(k-1)。其中G(p)表示fib数列模素数p的循环节长度，所以现在的问题是如何求fib数列模一个

小于10^6的素数p的循环节长度。

求fib数列模p(p是素数)的最小循环节方法：

暴力枚举fib[i]%p==0的最小的i，然后记录pos=i+1，设a为fib[i]%p==0的前一位数，即a=fib[i-1]

那么我们知道fib数列模p的最小循环节长度一定是pos*x，那么也就是说现在要求一个最小的数x，满足,

求出x后，那么问题就解决了，剩下的就是合并等等。

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
  return *p1++;
}

inline void read(int &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

inline ll Gcd(ll a,ll b){
  return b?Gcd(b,a%b):a;
}

inline ll Lcm(ll a,ll b){
  return a*b/Gcd(a,b);
}

inline ll Pow(ll a,ll b,ll p){
  ll ret=1;
  for (;b;b>>=1,a=a*a%p)
    if (b&1)
      ret=ret*a%p;
  return ret;
}

const int maxn=1000000;

int prime[maxn+5],num;
int vst[maxn+5];

inline void Pre(){
  for (int i=2;i<=maxn;i++){
    if (!vst[i]) prime[++num]=i;
    for (int j=1;j<=num && (ll)prime[j]*i<=maxn;j++){
      vst[prime[j]*i]=1;
      if (i%prime[j]==0)
	break;
    }
  }
}

inline ll calc(ll p){
  ll a=3,f1=1,f2=1,f3=2%p;  
  while (f3) f1=f2,f2=f3,f3=(f1+f2)%p,a++;
  ll ret=p-1;
  for (int i=1;(ll)i*i<=p-1;i++)
    if ((p-1)%i==0){
      if (Pow(f2,i,p)==1) ret=min(ret,(ll)i);
      if (Pow(f2,(p-1)/i,p)==1) ret=min(ret,(ll)(p-1)/i);
    }
  return ret*a;
}

inline ll Solve(ll n){
  int tem=n; ll ret=0;
  for (int i=1;prime[i]<=tem;i++)
    if (tem%prime[i]==0){
      ll tmp=1;
      while (tem%prime[i]==0) tem/=prime[i],tmp*=prime[i];
      tmp=tmp/prime[i]*calc(prime[i]);
      ret=!ret?tmp:Lcm(ret,tmp);
    }
  return ret;
}

int main(){
  int T,n,ca=0;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  read(T); Pre();
  while (T--){
    read(n);
    if (n==1)
      printf("Case #%d: %I64d\n",++ca,1LL);
    else
      printf("Case #%d: %I64d\n",++ca,Solve(n));
  }
  return 0;
}

法二：

来自：http://blog.csdn.net/acdreamers/article/details/10983813

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

    （1）把n素因子分解，即

    （2）分别计算Fib数模每个的循环节长度，假设长度分别是

    （3）那么Fib模n的循环节长度

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

 

     这里有一个优美的定理：Fib数模的最小循环节长度等于，其中表示Fib数模素数的最小循环节长度。可以看出我们现在最重要的就是求

 

 

对于求我们利用如下定理：

 

   如果5是模的二次剩余，那么循环节的的长度是的因子，否则，循环节的长度是的因子。

顺便说一句，对于小于等于5的素数，我们直接特殊判断，loop(2)=3,loop(3)=8,loop(5)=20。

那么我们可以先求出所有的因子，然后用矩阵快速幂来一个一个判断，这样时间复杂度不会很大。

 

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
  return *p1++;
}

inline void read(int &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

namespace F{
  struct Matrix{
    ll a,b,c,d;
    Matrix(ll a=0,ll b=0,ll c=0,ll d=0):a(a),b(b),c(c),d(d) { }
  };
  inline Matrix Mul(Matrix A,Matrix B,ll p){
    return Matrix((A.a*B.a+A.b*B.c)%p,(A.a*B.b+A.b*B.d)%p,(A.c*B.a+A.d*B.c)%p,(A.c*B.b+A.d*B.d)%p);  
  }
  inline Matrix Pow(Matrix a,ll b,ll p){
    Matrix ret(1,0,0,1);
    for (;b;b>>=1,a=Mul(a,a,p))
      if (b&1)
	ret=Mul(ret,a,p);
    return ret;
  }
  inline ll f(ll n,ll p){
    if (!n) return 0;
    return Pow(Matrix(1,1,1,0),n-1,p).a;
  }
}

inline ll Gcd(ll a,ll b){
  return b?Gcd(b,a%b):a;
}

inline ll Lcm(ll a,ll b){
  return a*b/Gcd(a,b);
}

inline ll Pow(ll a,ll b,ll p){
  ll ret=1;
  for (;b;b>>=1,a=a*a%p)
    if (b&1)
      ret=ret*a%p;
  return ret;
}

const int maxn=1000000;

int prime[maxn+5],num;
int vst[maxn+5];

inline void Pre(){
  for (int i=2;i<=maxn;i++){
    if (!vst[i]) prime[++num]=i;
    for (int j=1;j<=num && (ll)prime[j]*i<=maxn;j++){
      vst[prime[j]*i]=1;
      if (i%prime[j]==0)
	break;
    }
  }
}

int d[1005],cnt;  
int pri[25],q[25],len;  
  
void dfs(ll t,ll cur){  
  if(t==len+1) { d[++cnt]=cur; return; }    
  dfs(t+1,cur);  
  for (int i=1;i<=q[t];i++)  
    dfs(t+1,cur*=pri[t]);  
}  

inline void Mac(int n){  
  int tem=n; len=0;  
  for (int i=1;i<=num && (ll)prime[i]*prime[i]<=tem;i++)  
    if (tem%prime[i]==0){  
      pri[++len]=prime[i]; q[len]=0;  
      while (tem%prime[i]==0) tem/=prime[i],q[len]++;  
    }  
  if (tem!=1) pri[++len]=tem,q[len]=1;  
  cnt=0;  
  dfs(1,1);  
}

inline ll legendre(ll a,ll p){  
  if(Pow(a,(p-1)>>1,p)==1) return 1; return -1;  
}  

const ll ans[]={0,0,3,8,0,20};

inline ll calc(ll p){
  if (p<=5) return ans[p];
  ll n;
  if (legendre(5,p)==1) Mac(p-1); else Mac((p+1)<<1);
  sort(d+1,d+cnt+1);
  for (int i=1;i<=cnt;i++){
    if (F::f(d[i],p)==0 && F::f(d[i]+1,p)==1)
      return d[i];
  }
}

inline ll Solve(ll n){
  int tem=n; ll ret=0;
  for (int i=1;prime[i]<=tem;i++)
    if (tem%prime[i]==0){
      ll tmp=1;
      while (tem%prime[i]==0) tem/=prime[i],tmp*=prime[i];
      tmp=tmp/prime[i]*calc(prime[i]);
      ret=!ret?tmp:Lcm(ret,tmp);
    }
  return ret;
}



int main(){
  int T,n,ca=0;
  freopen("t.in","r",stdin);
  freopen("t2.out","w",stdout);
  read(T); Pre();
  while (T--){
    read(n);
    if (n==1)
      printf("Case #%d: %I64d\n",++ca,1LL);
    else
      printf("Case #%d: %I64d\n",++ca,Solve(n));
  }
  return 0;
}