ZOJ 3729 Arnold

1、这道题官方题解是用Baby-Step，Giant-Step，但是我写了三天都失败了，最后经ACdreamer指点，用“求斐波那契数模n的循环节”的方法才做出来的。

2、将题目中的坐标用数字代入，可以发现坐标组成了一个斐波那契数列。我们要求的就是斐波那契数列模N的循环节，求出这个循环节后，再除以二，就是本题的答案了。

3、代码是我参照ACdreamer博客上的标程，然后自己实现的，欢迎指正。

4、摘抄ACdreamer给出的方法，证明先不用管，需要用到Fermat定理及其推论，注意到代码中出现了勒让德符号（Legendre）：

对于一个正整数n，我们求Fib数模n的循环节的长度的方法如下：

    （1）把n素因子分解，即

    （2）分别计算Fib数模每个的循环节长度，假设长度分别是

    （3）那么Fib模n的循环节长度

从上面三个步骤看来，貌似最困难的是第二步，那么我们如何求Fib模的循环节长度呢？

 

     这里有一个优美的定理：Fib数模的最小循环节长度等于，其中表示Fib数模素数的最小循环节长度。可以看出我们现在最重要的就是求

 

 

对于求我们利用如下定理：

 

   如果5是模的二次剩余，那么循环节的的长度是的因子，否则，循环节的长度是的因子。

顺便说一句，对于小于等于5的素数，我们直接特殊判断，loop(2)=3,loop(3)=8,loop(5)=20。

那么我们可以先求出所有的因子，然后用矩阵快速幂来一个一个判断，这样时间复杂度不会很大。

 

 

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
typedef unsigned long long LL;
const int MAXN=2;
const int MAXM=2;
const int N=4000005,M=20005;//这两个数字偏大了，有兴趣的朋友可以调成刚刚好。
bool vis[N];
LL prime[M],pri[M],fac[M],num[M],number,MOD,n;
struct Matrix{ //矩阵类，实现乘法运算和初始化运算，注意MOD是个全局变量。一定要在乘法过程中取模而不是乘法结束后整体取模，因为本题数据可以很大，会出错。
  int n,m;
  LL a[MAXN][MAXM];
 void clear(){
  n=MAXN;m=MAXM;
  memset(a,0,sizeof(a));
 }
 Matrix operator * (const Matrix &b) const{
  Matrix temp;
  temp.clear();
  temp.n=n;temp.m=b.m;
  for(int i=0;i<n;i++)
   for(int j=0;j<b.m;j++)
    for(int k=0;k<m;k++){
     temp.a[i][j]+=a[i][k]*b.a[k][j];
  temp.a[i][j]%=MOD;
 }
  return temp;
 }
};
Matrix A;
Matrix pow_mod(Matrix A,LL i,LL mod){ //矩阵快速幂，用来验证我们假设的循环节是否正确。（如果正确，在经过循环节次变换后，会得到x=0，y=1）
    if(i==0){
    Matrix u;
    u.clear();
    u.n=MAXN;u.m=MAXM;
    u.a[0][0]=u.a[1][1]=1%mod;
    return u;
 }
   MOD=mod;
   Matrix temp=pow_mod(A,i>>1,mod);
   temp=temp*temp;
   if(i&1) temp=temp*A;
    return temp;
}
LL pow_mod(LL a,LL p,LL n){ //普通的快速幂。
    if(p==0) return 1;
 LL ans=pow_mod(a,p/2,n);
 ans=ans*ans%n;
 if(p%2==1) ans=ans*a%n;
 return ans;
}
int solve(LL n){ //求出n的所有因子，保存在pri里，因子的个数保存在num里，也就是将n分解质因数。
    int cnt=0;
 LL m=(LL)sqrt(n+0.5);
    for(int i=0;prime[i]<=m;i++){
    if(n%prime[i]==0){
        int temp=0;
     pri[cnt]=prime[i];
     while(n%prime[i]==0){
         temp++;
      n/=prime[i];
     }
     num[cnt]=temp;
     cnt++;
    }
 }
 if(n>1){
    pri[cnt]=n;
    num[cnt]=1;
    cnt++;
 }
 return cnt;
}
int work(LL n){ //求出所有可能的循环节，保存在fac里。
    int c=0;
 LL m=(LL)sqrt(n+0.5);
 for(LL i=1;i<=m;i++)
  if(n%i==0){
     if(i*i==n) fac[c++]=i;
     else{
        fac[c++]=i;
     fac[c++]=n/i;
     }
  }
 return c;
}
LL gcd(LL a,LL b){
   return b==0?a:gcd(b,a%b);
}
LL lcm(LL a,LL b){
   return a/gcd(a,b)*b;
}
LL legendre(LL a,LL p){ //勒让德符号，本来应该返回-1，但是我的LL是unsigned long long，故返回0。
   if(pow_mod(a,(p-1)/2,p)==1) return 1;
   else return 0;
}
LL find(LL n){
    int cnt=solve(n);
 LL ans=1;
 for(int i=0;i<cnt;i++){
  LL period=1;
  int c;
  if(pri[i]==2) period=3;
  else if(pri[i]==3) period=8;
  else if(pri[i]==5) period=20;
  else{
      if(legendre(5,pri[i])==1)
    c=work(pri[i]-1);
   else
    c=work(2*(pri[i]+1));
   sort(fac,fac+c);
   for(int k=0;k<c;k++){
    Matrix a=pow_mod(A,fac[k],pri[i]);
    LL x=(a.a[0][0]*0+a.a[0][1]*1)%pri[i];//用矩阵快速幂来加速求Fibonacci数列（mod N）
    LL y=(a.a[1][0]*0+a.a[1][1]*1)%pri[i]; 
    if(x==0&&y==1){
     period=fac[k];
     break;
    }
   }
  }
  for(LL k=1;k<num[i];k++)
   period*=pri[i];
  ans=lcm(ans,period);
 }
 return ans;//ans即为Fibonacci数列的循环节
}
void sieve(LL n){//筛法求素数 
    LL m=(LL)sqrt(n+0.5);
 memset(vis,0,sizeof(vis));
 for(LL i=2;i<=m;i++) if(!vis[i])
  for(LL j=i*i;j<=n;j+=i) vis[j]=true;
}
int generate_prime(LL n){ //得到素数表
    sieve(n);
    int c=0;
 for(LL i=2;i<=n;i++) if(!vis[i])
  prime[c++]=i;
 return c;
}
int main(){
 number=generate_prime(150000);
 A.clear();
 A.a[0][0]=0;
 A.a[0][1]=A.a[1][0]=A.a[1][1]=1;//A矩阵是在矩阵快速幂中用于加速求Fibonacci数列的。
 while(cin>>n){
  cout<<find(n)/2<<endl;//得出循环节后，除以二就是答案。至于为什么除以二，只要将题目的坐标用实际的数字代进去就能发现规律了。
 }
 return 0;
}