Fibonacc 数列模n的循环节

只谈解法,不说(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<iostream>#include<cstdio>#include<cmath>#include<cstring>#include<algorithm>#include<cstdlib>#include<vector>#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;k<time[i];k++)             rec*=num[i];        ans=lcm(ans,rec);    }        return ans;}int main() {    ll n;    init();    Euler_prime();    cin>>n;//  while(cin>>n) //  {        cout<<find_loop(n)<<endl;//  }    return 0;}