HDU 3978 斐波那契循环节

题意：给出f(f(f...f(n)...)) 总共嵌套k次。问最后模p的值是多少。

首先应该明白的是这个题有循环节的。一个数模N的循环节就是这个数分解成素因子乘积的形式p1^a1*p2^a2*p3^a3...后，斐波那契模pi^ai的循环节的最大公约数。

那么一个素数的k次幂的循环节=斐波那契模上这个素数的循环节乘上p^(k-1)。

而一个素数p的循环节 如果p>5并且是5的二次剩余，那么循环节就是(p-1)的因子，否则就是2*(p+1)的因子。所以2 3 5 的时候需要特判一下。

知道这些就能求每一次嵌套的循环节了，通过矩阵连乘即可得出答案。

#include <iostream>#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>using namespace std;#define maxn 20050bool isprime[maxn];int prime[maxn],nprime,mod[maxn],M;int primeloop[maxn];void getprime(){    long long i,j;    memset(isprime,1,sizeof(isprime));    nprime=0;    for(i=2; i<maxn; i++)        if(isprime[i])        {            prime[nprime++]=i;            for(j=i*i; j<maxn; j+=i) isprime[j]=0;        }}int factor[100][2],tol,fac[100][2],numfac;void findfac(int n,int f[][2],int &t){    int x=n,l=(int)sqrt(1.0*n);    t=0;    for(int i=0; i<50; i++) f[i][0]=f[i][1]=0;    for(int i=0; prime[i]<=l; i++)        if(x%prime[i]==0)        {            f[t][0]=prime[i];            while(x%prime[i]==0) f[t][1]++,x/=prime[i];            t++;        }    if(x>1) f[t][0]=x,f[t++][1]++;}const int MAX=2;typedef struct{    long long m[MAX][MAX];} Matrix;Matrix P,I;Matrix matrixmul(Matrix a,Matrix b) //矩阵乘法{    int i,j,k;    Matrix c;    for (i=0; i<MAX; i++)        for (j=0; j<MAX; j++)        {            c.m[i][j]=0;            for(k=0; k<MAX; k++)                c.m[i][j]+=((a.m[i][k]%M)*(b.m[k][j]%M))%M;            c.m[i][j]%=M;        }    return c;}void quickpow(int n,int &x,int &y){    Matrix m=P,b=I;    while(n)    {        if(n&1) b=matrixmul(b,m);        n>>=1,m=matrixmul(m,m);    }    x=b.m[0][0],y=b.m[1][0];}int gcd(int a,int b){    return b?gcd(b,a%b):a;}int exp_mod(int a,int b,int c){    int ans=1;    a%=c;    while(b)    {        if(b&1) ans=ans*a%c;        b>>=1,a=a*a%c;    }    return ans;}int minloop,ff[maxn],numff;void dfs(int num,int s=1){    if(num==numfac)    {        ff[numff++]=s;        return;    }    for(int i=0; i<=fac[num][1]; i++)        dfs(num+1,s),s*=fac[num][0];}int getPrimeLoop(int p){    if(p==2) return 3;    if(p==3) return 8;    if(p==5) return 20;    M=p;    if(exp_mod(5,(p-1)>>1,p)==1) p--;    else p=2*p+2;    findfac(p,fac,numfac);    minloop=1e9;    numff=0;    dfs(0,1);    sort(ff,ff+numff);    int x,y;    for(int i=1; i<numff; i++)    {        quickpow(ff[i]-1,x,y);        if(x==0&&y==1) return ff[i];    }}int getLoop(int p,int k){    int ret;    if(p>19583)        ret=getPrimeLoop(p);    else        ret=primeloop[p];    for(int i=0; i<k-1; i++) ret*=p;    return ret;}int getmod(int n){    findfac(n,factor,tol);    int ret=1,tem;    for(int i=0; i<tol; i++)        tem=getLoop(factor[i][0],factor[i][1]),ret=ret/gcd(ret,tem)*tem;    return ret;}int main(){    int t,ca=0,p,n,k,x,y;    getprime();    P.m[0][0]=P.m[0][1]=P.m[1][0]=1,P.m[1][1]=0;    I.m[0][0]=I.m[1][1]=1,I.m[1][0]=I.m[0][1]=0;    for(int i=0; i<2250; i++)        primeloop[prime[i]]=getPrimeLoop(prime[i]);    scanf("%d",&t);    while(t--)    {        scanf("%d%d%d",&n,&k,&p);        mod[0]=p;        for(int i=1; i<=k; i++) mod[i]=getmod(mod[i-1]);        for(int i=k; i>=0; i--)        {            if(i<k) n%=mod[i+1];            M=mod[i],quickpow(n,x,y),n=x;        }        printf("Case #%d: %d\n",++ca,n);    }    return 0;}