1. 程式人生 > >[斐波那契數列 求解二次剩餘&二次模方程 BSGS] COGS 2114 [CodeChef FN]斐波那契數

[斐波那契數列 求解二次剩餘&二次模方程 BSGS] COGS 2114 [CodeChef FN]斐波那契數

我們可以考慮斐波那契的通項公式 換元可以轉化成一個二次方程 根據求根公式求出後 再用BSGS求出指數

注意分奇偶考慮

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#include<cmath>
#define cl(x) memset(x,0,sizeof(x))
using namespace std;
typedef long long ll;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
  return *p1++;
}

inline void read(ll &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b; 
}


namespace BSGS{
  inline ll Pow(ll a,ll b,ll p){  
    ll ret=1;  
    for (;b;(a*=a)%=p,b>>=1)  
      if (b&1)  
	(ret*=a)%=p;  
    return ret;  
  }  
  struct hashmap{  
#define MOD 1000007  
    ll num;  
    ll key[1000005],value[1000005],next[1000005];  
    ll head[MOD];  
    inline void insert(ll k,ll v){  
      ll ad=k%MOD; key[++num]=k; value[num]=v; next[num]=head[ad]; head[ad]=num;  
    }  
    inline ll operator[](ll k){  
      ll ad=k%MOD;  
      for (int p=head[ad];p;p=next[p]) if (key[p]==k) return value[p];  
      return -1;  
    }  
    inline void clear(){  
      cl(head); num=0;  
    }  
  }M;
  ll lastA=-1,lastP=-1;
  inline ll Solve(ll A,ll B,ll P,ll phi){  
    A%=P;  
    if (!A && !B) return 1;  
    if (!A) return 1LL<<60;
    ll m=sqrt((double)P)+1;
    if (lastA!=A || lastP!=P){
      M.clear();
      ll tmp=1;  
      M.insert(1,0);  
      for (int i=1;i<m;i++){  
	(tmp*=A)%=P;  
	if (M[tmp]==-1) M.insert(tmp,i);  
      }
      lastA=A; lastP=P;
    }
    ll inv=1,T=Pow(A,phi-m,P);  
    for (int k=0;k<m;k++){  
      ll i=M[B*inv%P];  
      if (i!=-1) return k*m+i;  
      (inv*=T)%=P;  
    }  
    return 1LL<<60;  
  }  
}


namespace ROOT{
  ll P;
  inline ll Pow(ll a,ll b){
    ll ret=1;
    for (;b;b>>=1,a=a*a%P) if (b&1) ret=ret*a%P;
    return ret;
  }
  inline ll legendre(ll a){
    return Pow(a,(P-1)>>1);
  }
  struct abcd{
    ll a,b,w; //a+b*sqrt(w)
    abcd(ll a=0,ll b=0,ll w=0):a(a),b(b),w(w) { }
    friend abcd operator *(abcd A,abcd B){
      return abcd((A.a*B.a%P+A.b*B.b%P*A.w%P)%P,(A.a*B.b%P+A.b*B.a%P)%P,A.w);
    }
  };
  inline abcd Pow(abcd a,int b){
    abcd ret=abcd(1,0,a.w);
    for (;b;b>>=1,a=a*a) if (b&1) ret=ret*a;
    return ret;
  }
  inline ll calc(ll n){
    //if (P==2) return 1;
    n%=P;
    if (!n) return 0;
    if (legendre(n)==P-1) return -1;
    ll a,w;
    while (1){
      a=rand()%P;
      w=((a*a-n)%P+P)%P;
      if (legendre(w)==P-1) break;
    }
    return Pow(abcd(a,1,w),(P+1)>>1).a;
  }
  inline pair<ll,ll> Root(ll n,ll p){
    P=p;
    ll ans=calc(n);
    if (ans==-1) return make_pair(-1LL,-1LL);
    return make_pair(ans,P-ans);
  }
}

ll P,A,SQRT5,INV2;

inline ll Pre(ll p){
  P=p;
  INV2=(P+1)>>1;
  SQRT5=ROOT::Root(5,p).first;
  A=(1+SQRT5)*INV2%P;
}

inline ll Pow(ll a,ll b){
  ll ret=1;
  for (;b;b>>=1,a=a*a%P) if (b&1) ret=ret*a%P;
  return ret;
}

inline ll Inv(ll a){
  return Pow(a,P-2);
}


inline ll SolveEven(ll F){
  ll C=F*SQRT5%P;
  pair<ll,ll> tem=ROOT::Root(C*C+4,P);
  if (tem.first==-1) return 1LL<<60;
  ll x1=(C+tem.first)*INV2%P,x2=(C+tem.second)*INV2%P;
  ll ret=1LL<<60;
  tem=ROOT::Root(x1,P);
  if (tem.first!=-1)
    ret=min(ret,BSGS::Solve(A,tem.first,P,P-1)<<1);
  if (tem.second!=tem.first)
    ret=min(ret,BSGS::Solve(A,tem.second,P,P-1)<<1);
  if (x1!=x2){
    tem=ROOT::Root(x2,P);
    if (tem.first!=-1)
      ret=min(ret,BSGS::Solve(A,tem.first,P,P-1)<<1);
    if (tem.second!=tem.first)
      ret=min(ret,BSGS::Solve(A,tem.second,P,P-1)<<1);
  }
  return ret;
}


inline ll SolveOdd(ll F){
  ll C=F*SQRT5%P;
  pair<ll,ll> tem=ROOT::Root(C*C+P-4,P);
  if (tem.first==-1) return 1LL<<60;
  ll x1=(C+tem.first)*INV2%P,x2=(C+tem.second)*INV2%P;
  ll ret=1LL<<60;
  tem=ROOT::Root(x1*A,P);
  if (tem.first!=-1)
    ret=min(ret,(BSGS::Solve(A,tem.first,P,P-1)<<1)-1);
  if (tem.second!=tem.first)
    ret=min(ret,(BSGS::Solve(A,tem.second,P,P-1)<<1)-1);
  if (x1!=x2){
    tem=ROOT::Root(x2*A,P);
    if (tem.first!=-1)
      ret=min(ret,(BSGS::Solve(A,tem.first,P,P-1)<<1)-1);
    if (tem.second!=tem.first)
      ret=min(ret,(BSGS::Solve(A,tem.second,P,P-1)<<1)-1);
  }
  return ret;
}

int main(){
  ll T,f,p,tem,Ans;
  freopen("fn.in","r",stdin);
  freopen("fn.out","w",stdout);
  read(T);
  while (T--){
    read(f); read(p);
    Pre(p);
    Ans=1LL<<60;
    tem=SolveEven(f);
    Ans=min(Ans,tem);
    tem=SolveOdd(f);
    Ans=min(Ans,tem);
    if (Ans==1LL<<60)
      printf("-1\n");
    else
      printf("%lld\n",Ans);
  }
  return 0;
}