首页 > 代码库 > 斐波那契数列寻找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 循环节 模板
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。