广义Fibonacci数列找循环节 (二次剩余)
**问题:**给定
,满足
,求
的循
环节长度。
这里只说做法
我们先写出递推式的特征式子 x^2 =ax + b,整理得到 x^2-ax-b=0求出 delta = a^2+4b
对于质因数小于等于delta的部分,我们选择暴力求循环节。
暴力求循环节的代码如下:
1int get_loop( LL p) //暴力得到不大于13的素数的循环节。
2{
3 LL a,b,c;
4 a = 0 ;
5 b = 1 ;
6 for ( int i = 2; ; i++)
7 {
8 c = (a+3*b%p)%p; //此处为递推式
9 a = b;
10 b = c;
11 if (a==0&&b==1) return i-1;
12 }
13}
通常做法是先暴力求出来之后,写进一个表里。
通常不会有很多项。
对于质因数大于delta的部分,我们判断每个质因数是否是delta的二次剩余,如果是,枚举(prime-1)的因子,否则枚举2*(prime+1)的因子。
贴一个板子,需要注意如果是多个函数嵌套的形式,我们只需要从外向里,依次求循环节即可。
1/* ***********************************************
2Author :111qqz
3Created Time :Mon 31 Oct 2016 08:22:17 PM CST
4File Name :code/hdu/4291_2.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;
33struct Mat
34{
35 LL mat[2][2];
36 void clear()
37 {
38 ms(mat,0);
39 }
40}M,M1;
41Mat mul (Mat a,Mat b,LL mod)
42{
43 Mat c;
44 c.clear();
45 for ( int i = 0 ; i < 2 ; i++)
46 for ( int j = 0 ; j < 2 ; j++)
47 for ( int k = 0 ; k < 2 ; k++)
48 c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j]%mod)%mod;
49 return c;
50}
51Mat mat_ksm(Mat a,LL b,LL mod)
52{
53 Mat res;
54 res.clear();
55 for ( int i = 0 ; i < 2 ; i++) res.mat[i][i] = 1;
56 while (b>0)
57 {
58 if (b&1) res = mul(res,a,mod);
59 b = b >> 1LL;
60 a = mul(a,a,mod);
61 }
62 return res;
63}
64LL gcd(LL a,LL b)
65{
66 return b?gcd(b,a%b):a;
67}
68const int N = 1E6+7;
69bool prime[N];
70int p[N];
71void isprime() //一个普通的筛
72{
73 int cnt = 0 ;
74 ms(prime,true);
75 for ( int i = 2 ; i < N ; i++)
76 {
77 if (prime[i])
78 {
79 p[cnt++] = i ;
80 for ( int j = i+i ; j < N ; j+=i) prime[j] = false;
81 }
82 }
83}
84LL ksm( LL a,LL b,LL mod)
85{
86 LL res = 1;
87 while (b>0)
88 {
89 if (b&1) res = (res * a) % mod;
90 b = b >> 1LL;
91 a = a * a % mod;
92 }
93 return res;
94}
95LL legendre(LL a,LL p) //勒让德符号,判断二次剩余
96{
97 if (ksm(a,(p-1)>>1,p)==1) return 1;
98 return -1;
99}
100LL pri[N],num[N];//分解质因数的底数和指数。
101int cnt; //质因子的个数
102void solve(LL n,LL pri[],LL num[])
103{
104 cnt = 0 ;
105 for ( int i = 0 ; p[i] * p[i] <= n ; i++)
106 {
107 if (n%p[i]==0)
108 {
109 int Num = 0 ;
110 pri[cnt] = p[i];
111 while (n%p[i]==0)
112 {
113 Num++;
114 n/=p[i];
115 }
116 num[cnt] = Num;
117 cnt++;
118 }
119 }
120 if (n>1)
121 {
122 pri[cnt] = n;
123 num[cnt] = 1;
124 cnt++;
125 }
126}
127LL fac[N];
128int cnt2; //n的因子的个数
129void get_fac(LL n)//得到n的所有因子
130{
131 cnt2 = 0 ;
132 for (int i = 1 ; i*i <= n ; i++)
133 {
134 if (n%i==0)
135 {
136 if (i*i!=n)
137 {
138 fac[cnt2++] = i ;
139 fac[cnt2++] = n/i;
140 }
141 else fac[cnt2++] = i;
142 }
143 }
144}
145int get_loop( LL p) //暴力得到不大于13的素数的循环节。
146{
147 LL a,b,c;
148 a = 0 ;
149 b = 1 ;
150 for ( int i = 2; ; i++)
151 {
152 c = (a+3*b%p)%p;
153 a = b;
154 b = c;
155 if (a==0&&b==1) return i-1;
156 }
157}
158/*
159 2->3
160 3->2
161 5->12
162 7->16
163 11->8
164 13->52
165 */
166const LL LOOP[10]={3,2,12,16,8,52};
167LL ask_loop(int id)
168{
169 return LOOP[id];
170}
171LL find_loop(LL n)
172{
173 solve(n,pri,num);
174 LL ans = 1;
175 for ( int i = 0 ; i < cnt ; i++)
176 {
177 LL rec = 1;
178 if (pri[i]<=13) rec = ask_loop(i);
179 else
180 {
181 if (legendre(13,pri[i])==1) //13就是delta值
182 get_fac(pri[i]-1);
183 else get_fac((pri[i]+1)*(3-1)); //为什么可以假设循环节不大于2*(p+1)???
184 sort(fac,fac+cnt2);
185 for ( int j = 0 ; j < cnt2 ; j++) //挨个验证因子
186 {
187 Mat tmp = mat_ksm(M,fac[j],pri[i]); //下标从0开始,验证fac[j]为循环节,应该看fib[0]==fib[fac[j]]和fib[1]==fib[fac[j]+1]是否成立
188 tmp = mul(tmp,M1,pri[i]);
189 if (tmp.mat[0][0]==0&&tmp.mat[1][0]==1)
190 {
191 rec = fac[j];
192 break;
193 }
194 }
195
196 }
197 for ( int j = 1 ; j < num[i] ; j++)
198 rec *=pri[i];
199 ans = ans / gcd(ans,rec) * rec;
200 }
201 return ans;
202}
203void init()
204{
205 M.clear();
206 M.mat[0][1] = M.mat[1][0] = 1;
207 M.mat[1][1] = 3;
208 M1.clear();
209 M1.mat[1][0] = 1;
210}
211LL n;
212LL loop0 = 1E9+7;
213LL loop1,loop2;
214int main()
215{
216 #ifndef ONLINE_JUDGE
217 freopen("code/in.txt","r",stdin);
218 #endif
219 /*
220 printf("%d\n",get_loop(2));
221 printf("%d\n",get_loop(3));
222 printf("%d\n",get_loop(5));
223 printf("%d\n",get_loop(7));
224 printf("%d\n",get_loop(11));
225 printf("%d\n",get_loop(13));
226 */
227 init();
228 isprime();
229 while (~scanf("%lld\n",&n))
230 {
231 if (n==0||n==1)
232 {
233 printf("%lld\n",n);
234 continue;
235 }
236 LL loop1 = find_loop(loop0);
237 LL loop2 = find_loop(loop1);
238// printf("loop1:%lld loop2:%lld\n",loop1,loop2);
239 LL cur = n;
240 Mat ans = mat_ksm(M,cur-1,loop2);
241 ans = mul(ans,M1,loop2);
242 cur = ans.mat[1][0];
243 if (cur!=0&&cur!=1)
244 {
245 Mat ans = mat_ksm(M,cur-1,loop1);
246 ans = mul(ans,M1,loop1);
247 cur = ans.mat[1][0];
248 }
249 if (cur!=0&&cur!=1)
250 {
251 Mat ans = mat_ksm(M,cur-1,loop0);
252 ans = mul(ans,M1,loop0);
253 cur = ans.mat[1][0];
254 }
255 printf("%lld\n",cur);
256
257 }
258
259#ifndef ONLINE_JUDGE
260 fclose(stdin);
261#endif
262 return 0;
263}
再补一个,更为常用的,求斐波那契循环节的板子:
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}