广义Fibonacci数列找循环节 (二次剩余)

**问题:**给定 ,满足 ,求 的循

环节长度。

原理见广义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)的因子。

贴一个板子,需要注意如果是多个函数嵌套的形式,我们只需要从外向里,依次求循环节即可。

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

再补一个,更为常用的,求斐波那契循环节的板子:

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