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。

暴力解法：

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
 
 
 
/* ***********************************************
Author :111qqz
Created Time :Mon 31 Oct 2016 02:47:01 PM CST
File Name :code/hdu/4291.cpp
************************************************ */
 
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <cmath>
#include <cstdlib>
#include <ctime>
#define fst first
#define sec second
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
#define ms(a,x) memset(a,x,sizeof(a))
typedef long long LL;
#define pi pair < int ,int >
#define MP make_pair
 
using namespace std;
const double eps = 1E-8;
const int dx4[4]={1,0,0,-1};
const int dy4[4]={0,-1,1,0};
const int inf = 0x3f3f3f3f;
const LL loop0 = 1E9+7;
const LL loop1=222222224,loop2=183120; //暴力得到每一层的循环节。
/*
LL find_loop(LL mod)
{
    LL a,b,c;
    a = 0 ;
    b = 1 ;
    for ( int i = 2 ; ; i++)
    {
        c = (a+3*b%mod)% mod;
        a = b;
        b = c;
//        cout<<"a:"<<a<<" b:"<<b<<endl;
        if (a==0&&b==1)
        {
            return i-1; // i的时候得到第i项，判断循环的时候是判断g[i-1]==g[0]&&g[i]==g[1]，所以循环节长度是i-1;
        }
    }
}*/
LL n;
struct Mat
{
    LL mat[2][2];
    void clear()
    {
        ms(mat,0);
    }
    void out()
    {
        cout<<mat[0][0]<<" "<<mat[0][1]<<endl;
        cout<<mat[1][0]<<" "<<mat[1][1]<<endl;
        cout<<endl;
    }
}M,M1;
Mat mul(Mat a,Mat b,LL MOD)
{
    Mat c;
    c.clear();
    for ( int i = 0 ; i < 2 ; i++)
        for  ( int j = 0 ; j < 2 ; j++)
            for ( int k = 0 ; k < 2 ; k++)
                c.mat[i][j] = (c.mat[i][j] + a.mat[i][k]*b.mat[k][j]%MOD)%MOD;
    return c;
}
Mat power(Mat a,LL b,LL MOD)
{
    Mat res;
    res.clear();
    for ( int i = 0 ;i < 2 ; i++) res.mat[i][i] = 1;
    while (b>0)
    {
        if (b&1) res = mul(res,a,MOD);
        b = b>>1LL;
        a = mul(a,a,MOD);
    }
    return res;
}
void Mat_init()
{
    M.clear();
    M.mat[0][1] = 1;
    M.mat[1][0] = 1;
    M.mat[1][1] = 3;
    M1.clear();
    M1.mat[1][0] = 1;
}
int main()
{
        #ifndef  ONLINE_JUDGE
        freopen("code/in.txt","r",stdin);
  #endif
        /*
        loop1 = find_loop(loop0);
        printf("loop1:%lld\n",loop1);
        loop2 = find_loop(loop1);
        printf("loop2:%lld\n",loop2);
        */
        Mat_init();
        while (scanf("%lld\n",&n)!=EOF)
        {
            if (n==0)
            {
                puts("0");
                continue;
            }
            if (n==1)
            {
                puts("1");
                continue;
            }
 
            LL cur = n;
            Mat ans;
            ans = power(M,cur-1,loop2);
            ans = mul(ans,M1,loop2);
            cur = ans.mat[1][0];
            if (cur!=0&&cur!=1) //三次快速幂，每次都要记得判断初始项。
            {
 
                ans = power(M,cur-1,loop1);
                ans = mul(ans,M1,loop1);
                cur = ans.mat[1][0];
            }
            if (cur!=0&&cur!=1)
            {
                ans = power(M,cur-1,loop0);
                ans = mul(ans,M1,loop0);
                cur = ans.mat[1][0];
            }
            printf("%lld\n",cur);
        }
 
 
#ifndef ONLINE_JUDGE  
        fclose(stdin);
#endif
        return 0;
}
 
 
 
 

再来一个比较一般的做法：

参考递推式循环节

Acdreamer的博客_广义Fibonacci数列找循环节

摘重点：

今天早上起来后，看了下代码，为什么要判断5是不是p的模二次剩余呢，为什么是5呢

想了想，5对于斐波那契数列来讲，不就是x^2=x+1的delta么，那么这题的递推式是x^2=3*x+1，delta=3*3+4=13，然后我就把勒让德符号判断二次剩余那里改成13，然后对应的暴力出13及13以内的素数对应的循环节，交了一发，AC了

 

 

      所以综上所述：

是模的二次剩余时，枚举的因子

是模的二次非剩余时，枚举的因子

对于第二种非剩余的情况，理论上是枚举(p+1)(p-1)的因子，实际上常常只是枚举2*(p+1)的因子

对于c=5的情况，是有定理：

screenshot-from-2016-10-31-21-23-45

不过对于c不等于5的情况。。。该结论是否一定成立呢…感觉不是很好证，求指教。

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

 

作者： CrazyKK

ex-ACMer@hust，stackoverflow-engineer@sensetime

说点什么

您将是第一位评论人！

提醒
wpDiscuz