斐波那契数列寻找mod n 循环节模板

该代码来自ACdreamer
  1 #include <iostream>
  2 #include <string.h>
  3 #include <algorithm>
  4 #include <stdio.h>
  5 #include <math.h>
  6 
  7 using namespace std;
  8 typedef unsigned long long LL;
  9 
 10 const int M = 2;
 11 
 12 struct Matrix
 13 {
 14     LL m[M][M];
 15 };
 16 
 17 Matrix A;
 18 Matrix I = {1,0,0,1};
 19 
 20 Matrix multi(Matrix a,Matrix b,LL MOD)
 21 {
 22     Matrix c;
 23     for(int i=0; i<M; i++)
 24     {
 25         for(int j=0; j<M; j++)
 26         {
 27             c.m[i][j] = 0;
 28             for(int k=0; k<M; k++)
 29                 c.m[i][j] = (c.m[i][j]%MOD + (a.m[i][k]%MOD)*(b.m[k][j]%MOD)%MOD)%MOD;
 30             c.m[i][j] %= MOD;
 31         }
 32     }
 33     return c;
 34 }
 35 
 36 Matrix power(Matrix a,LL k,LL MOD)
 37 {
 38     Matrix ans = I,p = a;
 39     while(k)
 40     {
 41         if(k & 1)
 42         {
 43             ans = multi(ans,p,MOD);
 44             k--;
 45         }
 46         k >>= 1;
 47         p = multi(p,p,MOD);
 48     }
 49     return ans;
 50 }
 51 
 52 LL gcd(LL a,LL b)
 53 {
 54     return b? gcd(b,a%b):a;
 55 }
 56 
 57 const int N = 400005;
 58 const int NN = 5005;
 59 
 60 LL num[NN],pri[NN];
 61 LL fac[NN];
 62 int cnt,c;
 63 
 64 bool prime[N];
 65 int p[N];
 66 int k;
 67 
 68 void isprime()
 69 {
 70     k = 0;
 71     memset(prime,true,sizeof(prime));
 72     for(int i=2; i<N; i++)
 73     {
 74         if(prime[i])
 75         {
 76             p[k++] = i;
 77             for(int j=i+i; j<N; j+=i)
 78                 prime[j] = false;
 79         }
 80     }
 81 }
 82 
 83 LL quick_mod(LL a,LL b,LL m)
 84 {
 85     LL ans = 1;
 86     a %= m;
 87     while(b)
 88     {
 89         if(b & 1)
 90         {
 91             ans = ans * a % m;
 92             b--;
 93         }
 94         b >>= 1;
 95         a = a * a % m;
 96     }
 97     return ans;
 98 }
 99 
100 LL legendre(LL a,LL p)
101 {
102     if(quick_mod(a,(p-1)>>1,p)==1) return 1;
103     else                           return -1;
104 }
105 
106 void Solve(LL n,LL pri[],LL num[])
107 {
108     cnt = 0;
109     LL t = (LL)sqrt(1.0*n);
110     for(int i=0; p[i]<=t; i++)
111     {
112         if(n%p[i]==0)
113         {
114             int a = 0;
115             pri[cnt] = p[i];
116             while(n%p[i]==0)
117             {
118                 a++;
119                 n /= p[i];
120             }
121             num[cnt] = a;
122             cnt++;
123         }
124     }
125     if(n > 1)
126     {
127         pri[cnt] = n;
128         num[cnt] = 1;
129         cnt++;
130     }
131 }
132 
133 void Work(LL n)
134 {
135     c = 0;
136     LL t = (LL)sqrt(1.0*n);
137     for(int i=1; i<=t; i++)
138     {
139         if(n % i == 0)
140         {
141             if(i * i == n) fac[c++] = i;
142             else
143             {
144                 fac[c++] = i;
145                 fac[c++] = n / i;
146             }
147         }
148     }
149 }
150 
151 LL find_loop(LL n)
152 {
153     Solve(n,pri,num);
154     LL ans=1;
155     for(int i=0; i<cnt; i++)
156     {
157         LL record=1;
158         if(pri[i]==2)
159             record=3;
160         else if(pri[i]==3)
161             record=8;
162         else if(pri[i]==5)
163             record=20;
164         else
165         {
166             if(legendre(5,pri[i])==1)
167                 Work(pri[i]-1);
168             else
169                 Work(2*(pri[i]+1));
170             sort(fac,fac+c);
171             for(int k=0; k<c; k++)
172             {
173                 Matrix a = power(A,fac[k]-1,pri[i]);
174                 LL x = (a.m[0][0]%pri[i]+a.m[0][1]%pri[i])%pri[i];
175                 LL y = (a.m[1][0]%pri[i]+a.m[1][1]%pri[i])%pri[i];
176                 if(x==1 && y==0)
177                 {
178                     record = fac[k];
179                     break;
180                 }
181             }
182         }
183         for(int k=1; k<num[i]; k++)
184             record *= pri[i];
185         ans = ans/gcd(ans,record)*record;
186     }
187     return ans;
188 }
189 
190 void Init()
191 {
192     A.m[0][0] = 1;
193     A.m[0][1] = 1;
194     A.m[1][0] = 1;
195     A.m[1][1] = 0;
196 }
197 
198 int main()
199 {
200     LL n;
201     Init();
202     isprime();
203     while(cin>>n)
204         cout<<find_loop(n)<<endl;
205     return 0;
206 }
斐波那契数列寻找mod n 循环节 模板

斐波那契数列寻找mod n 循环节模板
