hdu3977 - fibonacci模p的周期

假设N(m)为模m的周期，则有

N(ab)=lcm(N(a),N(b))，其中gcd(a,b)=1; /*这里证明使用了fibo数的性质，

N(p^k)=N(p)*p^(k-1)，其中p为质数   /*论文里这个没有证明，但是验证了10^14内的都满足

对于质数p,其周期可以推算得出：

当p%5=1或4时，它是p-1的因子

当p%4=2或3时，它是2*(p-1)的因子

如果便只用验证因子是否是循环节就行。

求fibo数，可以使用矩阵乘法，

#include<stdio.h>#include<string.h>#include<math.h>#define maxp 1000000#define LL long longLL P;int pri[maxp+10];int b[maxp+10];int tot;struct data {    int p,t;    LL r;} q[500];int qt;LL fab[15];void getprime(){    memset(b,0,sizeof(b));    tot = 0;    int i=2;    while (i<=maxp) {        while (b[i] && i<=maxp)            i++;        pri[++tot]=i;        int j=i;        while (j<=maxp) {            b[j]=1;            j+=i;        }    }    tot--;    return ;}LL pow_mod(LL x,LL y,LL z){    LL res = 1;    x = x%z;    while (y) {        if (y&1)            res = (res*x)%z;        x = (x*x)%z;        y>>=1;    }    return res;}struct mat {    LL a[3][3];};mat mat_mul(mat x,mat y,int mod){    mat z;    for (int i=1; i<=2; i++)        for (int j=1; j<=2; j++) {            z.a[i][j]=0;            for (int k=1; k<=2; k++)                z.a[i][j] = (z.a[i][j]+(x.a[i][k]*y.a[k][j])%mod)%mod;        }    return z;}mat mat_pow(mat x,int n,int mod){    mat res;    res.a[1][1]=res.a[2][2]=1;    res.a[1][2]=res.a[2][1]=0;    while (n) {        if (n&1)            res = mat_mul(res,x,mod);        x = mat_mul(x,x,mod);        n>>=1;    }    return res;}int get_fibo(int n,int mod){    if (n<=14)        return fab[n]%mod;    mat x;    x.a[1][1]=1;    x.a[1][2]=1;    x.a[2][1]=1;    x.a[2][2]=0;    x = mat_pow(x,n-2,mod);    return (x.a[1][2]+x.a[1][1])%mod;}LL getans(LL pt,LL pr){    for (LL i=1; i<=pt; i++) {        if (pt % i==0) {            if (get_fibo(i,pr)==0 && get_fibo(i-1,pr)==1)                return i;        }    }}LL calc(LL pr){    if (pr == 2)        return 3;    if (pr==3)        return 8;    if (pr==5)        return 20;    if (pow_mod(5,pr/2,pr)==1)        return getans(pr-1,pr);    else        return getans(2*pr+2,pr);}void fenjie(LL n,data* a,int& cnt){    LL m = n;    int i=1;    cnt = 0;    while (i<=tot && m>1) {        if (m%pri[i]==0) {            cnt++;            a[cnt].p = pri[i];            a[cnt].t = 0;            while (m%pri[i]==0)                a[cnt].t++,m/=pri[i];        }        i++;    }    if (m>1) {        cnt++;        a[cnt].p = m;        a[cnt].t = 1;    }    for (int i=1; i<=cnt; i++) {        a[i].r = calc(a[i].p);    }    return ;}LL gcd(LL x,LL y){    return y==0?x:gcd(y,x%y);}LL lcm(LL x,LL y){    return x/gcd(x,y)*y;}LL power(LL x,LL y){    LL res = 1;    while (y) {        if (y&1)            res *= x;        x =x*x;        y>>=1;    }    return res;}LL solve(){    LL res=1;    fenjie(P,q,qt);    for (int i=1; i<=qt; i++) {        LL tmp = q[i].r*power(q[i].p,q[i].t-1);        res = lcm(res,tmp);    }    return res;}int main(){    fab[1]=fab[2]=1;    for (int i=3; i<=13; i++)        fab[i]=fab[i-1]+fab[i-2];    getprime();    int cas,cast = 0;    scanf("%d",&cas);    while (cas--) {        scanf("%I64d",&P);        LL ans = solve();        printf("Case #%d: %I64d\n",++cast,ans);    }    return 0;}