递推式求循环节的方法

来源：互联网发布：一冠淘宝店铺转让编辑：程序博客网时间：2024/06/05 16:50

大致说下递推式循环节的解决方案：

problem: f(n)=a*f(n-1)+b*f(n-2),求f(n)%p的循环节

solution:

1.对p进行质因数分解，p = p1^a1 * p2^a2 * p3^a3 ... * pn^an

2.分别求 p1^a1,p2^a2,...,pn^an的循环节，然后取最小公倍数

至于怎么求对 px^ax 的循环节，这里用到了二次剩余

2.1 p mod px^ax 的循环节 = G(px) * px^(ax-1) , G(px) 就是 p mod px 的最小循环节

2.2 对于递推式，我们可以得到特征根方程 x^2=a*x+b ,delta=a*a+4*b

2.3 对于G(px),如果delta是模px的二次剩余，G(px)是px-1的因子，否则G(px)是(px-1)*(px+1)的因子

暴力找到最小的就可以了

eg: acdream oj 1124 喵喵的遗憾：http://acdream.info/problem?pid=1124

代码如下：

#include <cstdio>#include <iostream>#include <cmath>#include <cstring>#include <algorithm>using namespace std;#define maxn 35000typedef long long LL;int prime[maxn];void Prime(){    memset(prime,0,sizeof(prime));    for(int i=2;i<maxn;i++){        if(!prime[i]) prime[++prime[0]]=i;        for(int j=1;j<=prime[0]&&prime[j]<maxn/i;j++){            prime[prime[j]*i]=1;            if(i%prime[j]==0){                break;            }        }    }}LL factor[100][2];int fatcnt;int get_factors(LL n){    fatcnt=0;    LL tmp=n;    for(int i=1;prime[i]<=tmp/prime[i];i++){        factor[fatcnt][1]=0;        if(tmp%prime[i]==0){            factor[fatcnt][0]=prime[i];            while(tmp%prime[i]==0){                tmp/=prime[i];                factor[fatcnt][1]++;            }            fatcnt++;        }    }    if(tmp!=1){        factor[fatcnt][0]=tmp;        factor[fatcnt][1]=1;        fatcnt++;    }    return fatcnt;}LL gcd(LL a,LL b){    if(b==0){        return a;    }    else{        return gcd(b,a%b);    }}LL lcm(LL a,LL b){    return a/gcd(a,b)*b;}struct Matrix{    LL m[2][2];}E,D;Matrix Multi(Matrix A,Matrix B,LL mod){    Matrix ans;    for(int i=0;i<2;i++){        for(int j=0;j<2;j++){            ans.m[i][j]=0;            for(int k=0;k<2;k++){                ans.m[i][j]+=(A.m[i][k]*B.m[k][j])%mod;                if(ans.m[i][j]>=mod){                    ans.m[i][j]-=mod;                }            }        }    }    return ans;}void init(){    memset(E.m,0,sizeof(E.m));    memset(D.m,0,sizeof(D.m));    D.m[0][0]=D.m[0][1]=D.m[1][0]=1;    for(int i=0;i<2;i++){        E.m[i][i]=1;    }    Prime();}Matrix Pow(Matrix A,LL e,LL mod){    Matrix ans=E;    while(e){        if(e&1){            ans=Multi(ans,A,mod);        }        A=Multi(A,A,mod);        e>>=1;    }    return ans;}LL Pow(LL a,LL b,LL mod){    LL ans=1;    while(b){        if(b&1){            ans=(ans*a)%mod;        }        a=(a*a)%mod;        b>>=1;    }    return ans;}int legendre(LL a,LL p){    if(Pow(a,(p-1)>>1,p)==1){        return 1;    }    else{        return -1;    }}int f0=1,f1=1;LL get_fib(LL n,LL mod){    if(mod==1) return 0;    return Pow(D,n,mod).m[0][0]%mod;}LL fac[maxn],GG[maxn];LL G(LL p){    if(p<maxn && GG[p]!=-1) return GG[p];    LL num;    if(legendre(5,p)==1){        num=p-1;    }    else{        num=2*(p+1);    }    int cnt=0;    for(LL i=1;i*i<=num;i++){        if(num%i==0){            fac[cnt++]=i;            if(i*i!=num){                fac[cnt++]=num/i;            }        }    }    sort(fac,fac+cnt);    LL ans;    for(int i=0;i<cnt;i++){        if(get_fib(fac[i],p)==f0&&get_fib(fac[i]+1,p)==f1){            ans=fac[i];            break;        }    }    if(p<maxn) GG[p]=ans;    return ans;}LL find_loop(LL n){    get_factors(n);    LL ans=1;    for(int i=0;i<fatcnt;i++)    {        LL record=1;        if(factor[i][0]==2) record=3;        else if(factor[i][0]==3) record=8;        else if(factor[i][0]==5) record=20;        else  record=G(factor[i][0]);        for(int j=1;j<factor[i][1];j++)            record*=factor[i][0];        ans=lcm(ans,record);    }    return ans;} int main(){    init();    memset(GG,-1,sizeof(GG));    int T;    LL N,P;    scanf("%d",&T);    while(T--)    {        scanf("%lld%lld",&N,&P);        LL mod1=P;        LL mod2=find_loop(mod1);        LL mod3=find_loop(mod2);        N=get_fib(N,mod3);        N=get_fib(N,mod2);        N=get_fib(N,mod1);        printf("%lld\n",N);    }    return 0;}

0 0