Fib数模n的循环节
我們知道Fibonacci數列,現在我們來求一個Fib數模n的循環節的長度。
對于一個正整數n,我們求Fib數模n的循環節的長度的方法如下:
? ? (1)把n素因子分解,即
? ? (2)分別計算Fib數模每個的循環節長度,假設長度分別是
? ? (3)那么Fib模n的循環節長度
從上面三個步驟看來,貌似最困難的是第二步,那么我們如何求Fib模的循環節長度呢?
?
?????這里有一個優美的定理:Fib數模的最小循環節長度等于,其中表示Fib數模素數的最小循環節長度。可以看出我們現在最重要的就是求
?
?
對于求我們利用如下定理:
?
???如果5是模的二次剩余,那么循環節的的長度是的因子,否則,循環節的長度是的因子。
順便說一句,對于小于等于5的素數,我們直接特殊判斷,loop(2)=3,loop(3)=8,loop(5)=20。
那么我們可以先求出所有的因子,然后用矩陣快速冪來一個一個判斷,這樣時間復雜度不會很大。
?
模板代碼:
#include <iostream> #include <string.h> #include <algorithm> #include <stdio.h> #include <math.h>using namespace std; typedef unsigned long long LL;const int M = 2;struct Matrix {LL m[M][M]; };Matrix A; Matrix I = {1,0,0,1};Matrix multi(Matrix a,Matrix b,LL MOD) {Matrix c;for(int i=0; i<M; i++){for(int j=0; j<M; j++){c.m[i][j] = 0;for(int k=0; k<M; k++)c.m[i][j] = (c.m[i][j]%MOD + (a.m[i][k]%MOD)*(b.m[k][j]%MOD)%MOD)%MOD;c.m[i][j] %= MOD;}}return c; }Matrix power(Matrix a,LL k,LL MOD) {Matrix ans = I,p = a;while(k){if(k & 1){ans = multi(ans,p,MOD);k--;}k >>= 1;p = multi(p,p,MOD);}return ans; }LL gcd(LL a,LL b) {return b? gcd(b,a%b):a; }const int N = 400005; const int NN = 5005;LL num[NN],pri[NN]; LL fac[NN]; int cnt,c;bool prime[N]; int p[N]; int k;void isprime() {k = 0;memset(prime,true,sizeof(prime));for(int i=2; i<N; i++){if(prime[i]){p[k++] = i;for(int j=i+i; j<N; j+=i)prime[j] = false;}} }LL quick_mod(LL a,LL b,LL m) {LL ans = 1;a %= m;while(b){if(b & 1){ans = ans * a % m;b--;}b >>= 1;a = a * a % m;}return ans; }LL legendre(LL a,LL p) {if(quick_mod(a,(p-1)>>1,p)==1) return 1;else return -1; }void Solve(LL n,LL pri[],LL num[]) {cnt = 0;LL t = (LL)sqrt(1.0*n);for(int i=0; p[i]<=t; i++){if(n%p[i]==0){int a = 0;pri[cnt] = p[i];while(n%p[i]==0){a++;n /= p[i];}num[cnt] = a;cnt++;}}if(n > 1){pri[cnt] = n;num[cnt] = 1;cnt++;} }void Work(LL n) {c = 0;LL t = (LL)sqrt(1.0*n);for(int i=1; i<=t; i++){if(n % i == 0){if(i * i == n) fac[c++] = i;else{fac[c++] = i;fac[c++] = n / i;}}} }LL find_loop(LL n) {Solve(n,pri,num);LL ans=1;for(int i=0; i<cnt; i++){LL record=1;if(pri[i]==2)record=3;else if(pri[i]==3)record=8;else if(pri[i]==5)record=20;else{if(legendre(5,pri[i])==1)Work(pri[i]-1);elseWork(2*(pri[i]+1));sort(fac,fac+c);for(int k=0; k<c; k++){Matrix a = power(A,fac[k]-1,pri[i]);LL x = (a.m[0][0]%pri[i]+a.m[0][1]%pri[i])%pri[i];LL y = (a.m[1][0]%pri[i]+a.m[1][1]%pri[i])%pri[i];if(x==1 && y==0){record = fac[k];break;}}}for(int k=1; k<num[i]; k++)record *= pri[i];ans = ans/gcd(ans,record)*record;}return ans; }void Init() {A.m[0][0] = 1;A.m[0][1] = 1;A.m[1][0] = 1;A.m[1][1] = 0; }int main() {LL n;Init();isprime();while(cin>>n)cout<<find_loop(n)<<endl;return 0; }
典型題目
題目一:http://acm.hdu.edu.cn/showproblem.php?pid=3977
題目二:http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=5075
題目三:http://acm.hdu.edu.cn/showproblem.php?pid=3978
題目四:http://acm.hdu.edu.cn/showproblem.php?pid=4291
上面的題目一和題目二基本上是模板題,沒有什么可以說的。對于第三題和第四題,我們可以看出題目四實際上是
題目三的簡單版,都是一層一層找循環節,我們找出每一層的循環節后即可解決。
下面給出題目三的代碼:
#include <iostream> #include <string.h> #include <algorithm> #include <stdio.h> #include <math.h>using namespace std; typedef long long LL;const int M=2;struct Matrix {LL m[M][M]; };Matrix per= {1,0,0,1};Matrix multi(Matrix a,Matrix b,LL MOD) {Matrix c;int i,j,k;for(i=0; i<M; i++){for(j=0; j<M; j++){c.m[i][j]=0;for(k=0; k<M; k++){c.m[i][j]+=a.m[i][k]*b.m[k][j]%MOD;}c.m[i][j]%=MOD;}}return c; }Matrix power(Matrix a,LL k,LL MOD) {Matrix ans=per,p=a;while(k){if(k&1){ans=multi(ans,p,MOD);k--;}k>>=1;p=multi(p,p,MOD);}return ans; }LL gcd(LL a,LL b) {return b? gcd(b,a%b):a; }LL quick_mod(LL a,LL b,LL m) {LL ans=1;a%=m;while(b){if(b&1){ans=ans*a%m;b--;}b>>=1;a=a*a%m;}return ans; }//勒讓德符號 int legendre(int a,int p) {if(quick_mod(a,(p-1)>>1,p)==1) return 1;else return -1; }const int N=1000005; const int NN=50005;bool prime[N]; int p[N]; int num[NN],pri[NN]; int num1[NN],pri1[NN]; int arr[NN]; int loop[N]; int k,cnt,c;void isprime() {k=0;int i,j;memset(prime,true,sizeof(prime));for(i=2; i<N; i++){if(prime[i]){p[k++]=i;for(j=i+i; j<N; j+=i){prime[j]=false;}}} }void find(int n,int pri[],int num[]) {cnt=0;int t=(int)sqrt(1.0*n);for(int i=0; p[i]<=t; i++){if(n%p[i]==0){int a=0;pri[cnt]=p[i];while(n%p[i]==0){a++;n/=p[i];}num[cnt]=a;cnt++;}}if(n>1){pri[cnt]=n;num[cnt]=1;cnt++;} }void dfs(int dept,int product=1) {if(dept==cnt){arr[c++]=product;return;}for(int i=0; i<=num1[dept]; i++){dfs(dept+1,product);product*=pri1[dept];} }int find_loop(int n) {find(n,pri,num);int cnt1=cnt;LL ans=1;for(int i=0; i<cnt1; i++){c=0;int record=1;if(pri[i]==2)record=3;else if(pri[i]==3)record=8;else if(pri[i]==5)record=20;else{if(legendre(5,pri[i])==1)find(pri[i]-1,pri1,num1);elsefind(2*(pri[i]+1),pri1,num1);dfs(0,1);sort(arr,arr+c);for(int k=0; k<c; k++){Matrix A;A.m[0][0]=1;A.m[0][1]=1;A.m[1][0]=1;A.m[1][1]=0;Matrix a=power(A,arr[k]-1,pri[i]);int x=(a.m[0][0]+a.m[0][1])%pri[i];int y=(a.m[1][0]+a.m[1][1])%pri[i];if(x==1&&y==0){record=arr[k];break;}}}for(int k=1; k<num[i]; k++)record*=pri[i];ans=ans/gcd(ans,record)*record;}return ans; }void Solve(int p,int k) {loop[0]=p;for(int i=1; i<=k; i++)loop[i]=find_loop(loop[i-1]); }int work(int n,int k,int p) {int t=n;LL ret,MOD;Matrix ans;Matrix A;A.m[0][0]=1;A.m[0][1]=1;A.m[1][0]=1;A.m[1][1]=0;Solve(p,k);for(int i=k; i>=0; i--){MOD=loop[i];ans=power(A,t,MOD);ret=(ans.m[1][0]+ans.m[1][1])%MOD;t=ret;}return ret; }int main() {isprime();int T,n,k,p,tt=1;scanf("%d",&T);while(T--){scanf("%d%d%d",&n,&k,&p);printf("Case #%d: %d\n",tt++,work(n,k,p));}return 0; }
?
題目:http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1195
?
分析:由于本題數據特別多,開始一直TLE,后來發現矩陣乘法里面如果取模太多會導致速度變得很慢,以前POJ上
???? 的一道題關于矩陣求和也是因為這個TLE,盡量減少一些取模,速度會有幾十倍甚至幾百倍的提升。
?
代碼:
#include <iostream> #include <string.h> #include <algorithm> #include <stdio.h> #include <math.h> #include <map>using namespace std; typedef long long LL;const int M = 2;struct Matrix {LL m[M][M]; };Matrix A; Matrix I = {1,0,0,1};Matrix multi(Matrix a,Matrix b,LL MOD) {Matrix c;for(int i=0; i<M; i++){for(int j=0; j<M; j++){c.m[i][j] = 0;for(int k=0; k<M; k++)c.m[i][j] += a.m[i][k] * b.m[k][j];c.m[i][j] %= MOD;}}return c; }Matrix power(Matrix a,LL k,LL MOD) {Matrix ans = I,p = a;while(k){if(k & 1){ans = multi(ans,p,MOD);k--;}k >>= 1;p = multi(p,p,MOD);}return ans; }int gcd(int a,int b) {return b? gcd(b,a%b):a; }const int N = 400005; const int NN = 5005;int num[NN],pri[NN]; int fac[NN]; bool flag[NN]; int c;bool prime[N]; int p[N]; int k;int cnt1; int pri1[NN],num1[NN];void isprime() {k = 0;memset(prime,true,sizeof(prime));for(int i=2; i<N; i++){if(prime[i]){p[k++] = i;for(int j=i+i; j<N; j+=i)prime[j] = false;}} }LL quick_mod(LL a,LL b,LL m) {LL ans = 1;a %= m;while(b){if(b & 1){ans = ans * a % m;b--;}b >>= 1;a = a * a % m;}return ans; }int legendre(int a,int p) {if(quick_mod(a,(p-1)>>1,p)==1) return 1;else return -1; }void Solve(int n,int pri[],int num[],int &cnt) {cnt = 0;int t = (int)sqrt(1.0*n);for(int i=0; p[i]<=t; i++){if(n%p[i]==0){int a = 0;pri[cnt] = p[i];while(n%p[i]==0){a++;n /= p[i];}num[cnt] = a;cnt++;}}if(n > 1){pri[cnt] = n;num[cnt] = 1;cnt++;} }void dfs(int dept,int cnt,LL product,int pri1[],int num1[]) {if(dept == cnt){fac[c++] = product;return;}for(int i=0; i<=num1[dept]; i++){dfs(dept+1,cnt,product,pri1,num1);product *= pri1[dept];} }map<int,int> mp;LL find_loop(LL n) {int cnt = 0;Solve(n,pri,num,cnt);LL ans = 1;for(int i=0; i<cnt; i++){int record=1;if(mp.find(pri[i]) != mp.end()){record = mp[pri[i]];goto Test;}if(pri[i]==2)record=3;else if(pri[i]==3)record=8;else if(pri[i]==5)record=20;else{if(legendre(5,pri[i])==1){c = 0;Solve(pri[i]-1,pri1,num1,cnt1);dfs(0,cnt1,1,pri1,num1);}else{c = 0;Solve(2*(pri[i]+1),pri1,num1,cnt1);dfs(0,cnt1,1,pri1,num1);}sort(fac,fac+c);for(int r=0; r<c; r++)flag[r] = 1;for(int k=c-1; k >= 0; k--){if(!flag[k]) continue;Matrix a = power(A,fac[k]-1,pri[i]);int x = (a.m[0][0]%pri[i]+a.m[0][1]%pri[i])%pri[i];int y = (a.m[1][0]%pri[i]+a.m[1][1]%pri[i])%pri[i];if(x==1 && y==0){record = fac[k];}else{for(int j=0; j<=k; j++){if(fac[k] % fac[j] == 0)flag[j] = 0;}}}mp[pri[i]] = record;} Test:for(int k=1; k<num[i]; k++)record *= pri[i];ans=ans/gcd(ans,record)*record;}return ans; }void Init() {A.m[0][0] = 1;A.m[0][1] = 1;A.m[1][0] = 1;A.m[1][1] = 0; }int main() {int T,n;Init();mp.clear();isprime();scanf("%d",&T);while(T--){scanf("%d",&n);LL ans = find_loop(n);printf("%I64d\n",ans);}return 0; }
?
總結
以上是生活随笔為你收集整理的Fib数模n的循环节的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: HDU3208(区间指数和)
- 下一篇: HDU4405(概率DP求期望)