hdu3977 - fibonacci模p的周期

2019-04-14 16:08发布

站内文章 / 模拟电子
9310 0 1517

假设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 #include #include #define maxp 1000000 #define LL long long LL 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 ",++cast,ans); } return 0; }

热门文章