hdu 3977 Evil teacher (斐波那契数列循环节)

题目链接

题意:f[0] = 1,f[1] = 1,f[i] = f[i-1] + f[i-2] (i>=2),问最小的m满足f[n]%p==f[n+m]%p

思路:求斐波那契数列循环节。

参考了Acdreamer的博客_Fib数模n的循环节

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

(1)把n素因子分解,即

(2)分别计算Fib数模每个 的循环节长度,假设长度分别是

(3)那么Fib模n的循环节长度

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

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

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

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

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

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

这道题是模板题。博客中的模板代码很好理解…

换成了自己比较熟悉的矩阵构造方式,以及代码风格。

需要注意的是下标是从0开始的,当验证k是否为循环节的时候,应该验证f[k]和f[k+1]

  1/* ***********************************************
  2Author :111qqz
  3Created Time :Mon 31 Oct 2016 03:54:15 AM CST
  4File Name :code/hdu/3977.cpp
  5************************************************ */
  6
  7#include <cstdio>
  8#include <cstring>
  9#include <iostream>
 10#include <algorithm>
 11#include <vector>
 12#include <queue>
 13#include <set>
 14#include <map>
 15#include <string>
 16#include <cmath>
 17#include <cstdlib>
 18#include <ctime>
 19#define fst first
 20#define sec second
 21#define lson l,m,rt<<1
 22#define rson m+1,r,rt<<1|1
 23#define ms(a,x) memset(a,x,sizeof(a))
 24typedef long long LL;
 25#define pi pair < int ,int >
 26#define MP make_pair
 27
 28using namespace std;
 29const double eps = 1E-8;
 30const int dx4[4]={1,0,0,-1};
 31const int dy4[4]={0,-1,1,0};
 32const int inf = 0x3f3f3f3f;
 33LL P;
 34struct Mat
 35{
 36    LL mat[2][2];
 37    void clear()
 38    {
 39	ms(mat,0);
 40    }
 41}M,M1;
 42Mat mul (Mat a,Mat b,LL mod)
 43{
 44    Mat c;
 45    c.clear();
 46    for ( int i = 0 ; i < 2 ; i++)
 47	for ( int j = 0 ; j < 2 ; j++)
 48	    for ( int k  = 0 ; k < 2 ; k++)
 49		c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]%mod)%mod;
 50    return c;
 51}
 52Mat mat_ksm(Mat a,LL b,LL mod)
 53{
 54    Mat res;
 55    res.clear();
 56    for ( int i = 0 ; i < 2 ; i++) res.mat[i][i] = 1;
 57    while (b>0)
 58    {
 59	if (b&1) res = mul(res,a,mod);
 60	b = b >> 1LL;
 61	a = mul(a,a,mod);
 62    }
 63    return res;
 64}
 65LL gcd(LL a,LL b)
 66{
 67    return b?gcd(b,a%b):a;
 68}
 69const int N = 1E6+7;
 70bool prime[N];
 71int p[N];
 72void isprime() //一个普通的筛
 73{
 74    int cnt = 0 ;
 75    ms(prime,true);
 76    for ( int i = 2 ; i < N ; i++)
 77    {
 78	if (prime[i])
 79	{
 80	    p[cnt++] = i ;
 81	    for ( int j = i+i ; j < N ; j+=i) prime[j] = false;
 82	}
 83    }
 84}
 85LL ksm( LL a,LL b,LL mod)
 86{
 87   LL res = 1;
 88   while (b>0)
 89   {
 90       if (b&1) res = (res * a) % mod;
 91       b = b >> 1LL;
 92       a = a * a % mod;
 93   }
 94   return res;
 95}
 96LL legendre(LL a,LL p) //勒让德符号,判断二次剩余
 97{
 98    if (ksm(a,(p-1)>>1,p)==1) return 1;
 99    return -1;
100}
101LL pri[N],num[N];//分解质因数的底数和指数。
102int cnt; //质因子的个数
103void solve(LL n,LL pri[],LL num[])
104{
105    cnt = 0 ;
106    for ( int  i = 0 ; p[i] * p[i] <= n ; i++)
107    {
108	if (n%p[i]==0)
109	{
110	    int Num = 0 ;
111	    pri[cnt] = p[i];
112	    while (n%p[i]==0)
113	    {
114		Num++;
115		n/=p[i];
116	    }
117	    num[cnt] = Num;
118	    cnt++;
119	}
120    }
121    if (n>1)
122    {
123	pri[cnt] = n;
124	num[cnt] = 1;
125	cnt++;
126    }
127}
128LL fac[N];
129int cnt2; //n的因子的个数
130void get_fac(LL n)//得到n的所有因子
131{
132    cnt2 = 0 ;
133    for (int i =  1 ; i*i <= n ; i++)
134    {
135	if (n%i==0)
136	{
137	    if (i*i!=n)
138	    {
139		fac[cnt2++] = i ;
140		fac[cnt2++] = n/i;
141	    }
142	    else fac[cnt2++] = i;
143	}
144    }
145}
146LL find_loop(LL n)
147{
148    solve(n,pri,num);
149    LL ans = 1;
150    for ( int i = 0 ; i < cnt ; i++)
151    {
152	LL rec = 1;
153	if (pri[i]==2) rec = 3;
154	else if (pri[i]==3) rec = 8;
155	else if (pri[i]==5) rec = 20;
156	else
157	{
158	    if (legendre(5,pri[i])==1)
159		get_fac(pri[i]-1);
160	    else get_fac(2*pri[i]+2);
161	    sort(fac,fac+cnt2);
162	    for ( int j = 0 ; j < cnt2 ; j++) //挨个验证因子
163	    {
164		Mat tmp = mat_ksm(M,fac[j],pri[i]); //下标从0开始,验证fac[j]为循环节,应该看fib[0]==fib[fac[j]]和fib[1]==fib[fac[j]+1]是否成立
165		tmp = mul(tmp,M1,pri[i]);
166		if (tmp.mat[0][0]==1&&tmp.mat[1][0]==1)
167		{
168		    rec = fac[j];
169		    break;
170		}
171	    }
172
173	}
174	for ( int j = 1 ; j < num[i] ; j++)
175	    rec *=pri[i];
176	ans = ans / gcd(ans,rec) * rec;
177    }
178    return ans;
179}
180void init()
181{
182    M.clear();
183    M.mat[0][1] = M.mat[1][0] = M.mat[1][1] = 1;
184    M1.clear();
185    M1.mat[0][0] = M1.mat[1][0] = 1;
186}
187int main()
188{
189	#ifndef  ONLINE_JUDGE 
190	freopen("code/in.txt","r",stdin);
191  #endif
192	int T;
193	int cas = 0 ;
194	isprime();
195	scanf("%d",&T);
196	while (T--)
197	{
198	    init();
199	    scanf("%lld",&P);
200	    printf("Case #%d: ",++cas);
201	    LL ret = find_loop(P);
202	    printf("%lld\n",ret);
203	}
204
205  #ifndef ONLINE_JUDGE  
206  fclose(stdin);
207  #endif
208    return 0;
209}