Fibonacc 数列模n的循环节

2019-04-14 12:47发布

站内文章 / 模拟电子
14807 0 1156
只谈解法,不说(hui)证明
步骤:
对于一个斐波拉契数列求模n的循环节:
1:任意一个n都可以分解为 priem[i]^time[i]
2:对于每一个prime[i]求出对它的循环节 g[i],则任意prime[i]^time[i]循环节的大小为:
k[i]= g[i]^(time[i]-1)
3:则对于模n的循环节为lcm(k[i]); 求解g[i]
二次剩余:对于一元二次同余方程有无的情况的讨论
判定 对于a%p 来说
if(pow(a,(p-1)/2)mod p == 1 return 1//则a是p的二次剩余 else return -1//不是 —— 称为 legendre符号
ps:二次剩余只对奇素数成立,对2需要特殊判定 而对于prime[i]来说, 如果5是prime[i]的二次剩余则循环节为 priem[i]-1中的因子其中一个
否则 为 (prime[i]+1)*2 中因子的一个 用矩阵乘法判断即可:当 x==0&&y==1 时为循环节(x,y含义看代码)
Code: #include #include #include #include #include #include #include #define fo(i,a,b) for(int i=a;i<=b;i++) #define fod(i,a,b) for(int i=a;i>=b;i--) using namespace std; typedef long long ll; const int N = 400005; struct Marix { int val[2][2]; }; Marix mutli(Marix a,Marix b,ll mod) { Marix c; for(int i=0;i<2;i++) for(int j=0;j<2;j++) { c.val[i][j]=0; for(int k=0;k<2;k++) { c.val[i][j]+=(a.val[i][k]*b.val[k][j]); c.val[i][j]%=mod; } } return c; } ll gcd(ll a,ll b){return b==0?a:gcd(b,a%b);} ll lcm(ll a,ll b){return a*b/gcd(a,b);} ll prime[N],tot=0; void Euler_prime() { bool check[N]; memset(check,0,sizeof(check)); for(int i=2;i<=N;i++) { if(!check[i]) prime[++tot]=i; for(int j=1;j<=tot;j++) { if(prime[j]*i>N) break; check[prime[j]*i]=1; if(i%prime[j]==0) break; } } } ll num[N],time[N],cnt=0; void solve(ll n) { memset(time,0,sizeof(time)); ll t=ll(sqrt(n*1.00)); for(int i=1;prime[i]<=t;i++) { if(n%prime[i]==0) { num[++cnt]=prime[i]; while(n%prime[i]==0) { time[cnt]++; n/=prime[i]; } } } if(n>1) { num[++cnt]=n; time[cnt]++; } } ll quick_pow(ll a,ll b,ll c) { ll ans=1; for(;b;b>>=1) { if(b&1) { ans*=a; ans%=c;} a*=a; a%=c; } return ans; } Marix I,M; void init() { M.val[0][0]=M.val[0][1]=M.val[1][0]=1; M.val[1][1]=0; I.val[0][0]=I.val[1][1]=1; I.val[0][1]=I.val[1][0]=0; } Marix quick_Mpow(Marix a,ll k,ll mod) { Marix c=I,p=a; for(;k;k>>=1) { if(k&1) c=mutli(c,p,mod); p=mutli(p,p,mod); } return c; } ll legendre(ll a,ll p) { if(quick_pow(a,(p-1)>>1,p)==1) return 1; else return -1; } ll fac[N],c; void find_fac(ll n) { c=0; ll t=(ll)(sqrt(n*1.0)); 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); ll ans=1; for(int i=1;i<=cnt;i++) { ll rec=1; if(num[i]==2) rec=3; else if(num[i]==3) rec=8; else if(num[i]==5) rec=20; else { if(legendre(5,num[i])==1) find_fac(num[i]-1); else find_fac(2*(num[i]+1)); sort(fac,fac+c); for(int k=1;k<=c;k++) { Marix a=quick_Mpow(M,fac[k]-1,num[i]); ll x=(a.val[0][0]%num[i]+a.val[0][1]%num[i])%num[i]; ll y=(a.val[1][0]%num[i]+a.val[1][1]%num[i])%num[i]; if(x==1&&y==0) { rec=fac[k]; break; } } } for(int k=1;kreturn ans; } int main() { ll n; init(); Euler_prime(); cin>>n; // while(cin>>n) // { cout<// } return 0; }

热门文章