首页 > 代码库 > 斐波那契数列寻找mod n 循环节 模板

斐波那契数列寻找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             else143             {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         else165         {166             if(legendre(5,pri[i])==1)167                 Work(pri[i]-1);168             else169                 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 循环节 模板