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

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

  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}