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}