hdu 4291 A Short problem (矩阵快速幂+广义斐波那契循环节||暴力找循环节)
题意：
Given n (1 <= n <= 1018), You should solve for
g(g(g(n))) mod 109 + 7 where
g(n) = 3g(n - 1) + g(n - 2)
g(1) = 1
g(0) = 0思路：找循环节。首先由于模数固定，可以暴力一下找到循环节。
得到1E9+7的循环节是222222224,222222224的循环节是183120.
然后三次矩阵快速幂就行了。
需要注意每次都要判断那一层的n是否为0和1。
暴力解法：
1/* ***********************************************
2Author :111qqz
3Created Time :Mon 31 Oct 2016 02:47:01 PM CST
4File Name :code/hdu/4291.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;
33const LL loop0 = 1E9+7;
34const LL loop1=222222224,loop2=183120; //暴力得到每一层的循环节。
35/*
36LL find_loop(LL mod)
37{
38 LL a,b,c;
39 a = 0 ;
40 b = 1 ;
41 for ( int i = 2 ; ; i++)
42 {
43 c = (a+3*b%mod)% mod;
44 a = b;
45 b = c;
46// cout<<"a:"<<a<<" b:"<<b<<endl;
47 if (a==0&&b==1)
48 {
49 return i-1; // i的时候得到第i项，判断循环的时候是判断g[i-1]==g[0]&&g[i]==g[1]，所以循环节长度是i-1;
50 }
51 }
52}*/
53LL n;
54struct Mat
55{
56 LL mat[2][2];
57 void clear()
58 {
59 ms(mat,0);
60 }
61 void out()
62 {
63 cout<<mat[0][0]<<" "<<mat[0][1]<<endl;
64 cout<<mat[1][0]<<" "<<mat[1][1]<<endl;
65 cout<<endl;
66 }
67}M,M1;
68Mat mul(Mat a,Mat b,LL MOD)
69{
70 Mat c;
71 c.clear();
72 for ( int i = 0 ; i < 2 ; i++)
73 for ( int j = 0 ; j < 2 ; j++)
74 for ( int k = 0 ; k < 2 ; k++)
75 c.mat[i][j] = (c.mat[i][j] + a.mat[i][k]*b.mat[k][j]%MOD)%MOD;
76 return c;
77}
78Mat power(Mat a,LL b,LL MOD)
79{
80 Mat res;
81 res.clear();
82 for ( int i = 0 ;i < 2 ; i++) res.mat[i][i] = 1;
83 while (b>0)
84 {
85 if (b&1) res = mul(res,a,MOD);
86 b = b>>1LL;
87 a = mul(a,a,MOD);
88 }
89 return res;
90}
91void Mat_init()
92{
93 M.clear();
94 M.mat[0][1] = 1;
95 M.mat[1][0] = 1;
96 M.mat[1][1] = 3;
97 M1.clear();
98 M1.mat[1][0] = 1;
99}
100int main()
101{
102 #ifndef ONLINE_JUDGE
103 freopen("code/in.txt","r",stdin);
104 #endif
105 /*
106 loop1 = find_loop(loop0);
107 printf("loop1:%lld\n",loop1);
108 loop2 = find_loop(loop1);
109 printf("loop2:%lld\n",loop2);
110 */
111 Mat_init();
112 while (scanf("%lld\n",&n)!=EOF)
113 {
114 if (n==0)
115 {
116 puts("0");
117 continue;
118 }
119 if (n==1)
120 {
121 puts("1");
122 continue;
123 }
124
125 LL cur = n;
126 Mat ans;
127 ans = power(M,cur-1,loop2);
128 ans = mul(ans,M1,loop2);
129 cur = ans.mat[1][0];
130 if (cur!=0&&cur!=1) //三次快速幂，每次都要记得判断初始项。
131 {
132
133 ans = power(M,cur-1,loop1);
134 ans = mul(ans,M1,loop1);
135 cur = ans.mat[1][0];
136 }
137 if (cur!=0&&cur!=1)
138 {
139 ans = power(M,cur-1,loop0);
140 ans = mul(ans,M1,loop0);
141 cur = ans.mat[1][0];
142 }
143 printf("%lld\n",cur);
144 }
145
146
147#ifndef ONLINE_JUDGE
148 fclose(stdin);
149#endif
150 return 0;
151}
再来一个比较一般的做法：
参考递推式循环节
Acdreamer的博客_广义Fibonacci数列找循环节
摘重点：
今天早上起来后，看了下代码，为什么要判断5是不是p的模二次剩余呢，为什么是5呢想了想，5对于斐波那契数列来讲，不就是x^2=x+1的delta么，那么这题的递推式是x^2=3x+1，delta=33+4=13，然后我就把勒让德符号判断二次剩余那里改成13，然后对应的暴力出13及13以内的素数对应的循环节，交了一发，AC了
所以综上所述：
是模
的二次剩余时，枚举
的因子
的因子
对于第二种非剩余的情况，理论上是枚举(p+1)(p-1)的因子，实际上常常只是枚举2(p+1)的因子*
对于c=5的情况，是有定理：
不过对于c不等于5的情况。。。该结论是否一定成立呢…感觉不是很好证，求指教。
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)
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}