BZOJ5104 : Fib数列

斐波那契数列满足$f(n-1)f(n+1)-f(n)^2=(-1)^n$。

枚举$-1$的符号，根据二次剩余即可求出最多$4$个可能的$f(n+1)$的值。

那么根据$f(n)$和$f(n+1)$，对矩阵做BSGS求出最小的$n$即可。

时间复杂度$O(sqrt{P}log P)$。

#include<cstdio>
#include<algorithm>
#include<map>
using namespace std;
typedef long long ll;
typedef pair<ll,ll>PI;
const int P=1000000009,K=32000;
ll X,ans=-1,g[K+5][2][2],w[K+5][2];map<PI,int>T;
inline ll po(ll a,ll b,ll P){ll t=1;for(;b;b>>=1,a=a*a%P)if(b&1)t=t*a%P;return t;}
ll ToneLLi_Shanks(ll n,ll p){
  if(!n)return 0;
  if(p==2)return n&1?1:-1;
  if(po(n,p>>1,p)!=1)return -1;
  if(p&2)return po(n,p+1>>2,p);
  int s=__builtin_ctzll(p^1);
  ll q=p>>s,z=2;
  for(;po(z,p>>1,p)==1;z++);
  ll c=po(z,q,p);
  ll r=po(n,q+1>>1,p);
  ll t=po(n,q,p),tmp;
  for(int m=s,i;t!=1;){
    for(i=0,tmp=t;tmp!=1;i++)tmp=tmp*tmp%p;
    for(;i<--m;)c=c*c%p;
    r=r*c%p;
    c=c*c%p;
    t=t*c%p;
  }
  return r;
}
inline void mul(ll a[][2],ll b[][2],ll f[][2]){
  static ll c[2][2];
  c[0][0]=(a[0][0]*b[0][0]+a[0][1]*b[1][0])%P;
  c[0][1]=(a[0][0]*b[0][1]+a[0][1]*b[1][1])%P;
  c[1][0]=(a[1][0]*b[0][0]+a[1][1]*b[1][0])%P;
  c[1][1]=(a[1][0]*b[0][1]+a[1][1]*b[1][1])%P;
  f[0][0]=c[0][0];
  f[0][1]=c[0][1];
  f[1][0]=c[1][0];
  f[1][1]=c[1][1];
}
void work(ll X,ll Y){
  T.clear();
  for(int i=0;i<K;i++){
    ll a=(g[i][0][0]*X+g[i][0][1]*Y)%P;
    ll b=(g[i][1][0]*X+g[i][1][1]*Y)%P;
    T[PI(a,b)]=i;
  }
  for(int i=1;i<=K;i++){
    map<PI,int>::iterator it=T.find(PI(w[i][0],w[i][1]));
    if(it!=T.end()){
      ll now=i*K-it->second;
      if(ans==-1||ans>now)ans=now;
      return;
    }
  }
}
void solve(ll A,ll B,ll C){
  ll d=((B*B-4*A*C)%P+P)%P;
  d=ToneLLi_Shanks(d,P);
  if(d<0)return;
  d=(d%P+P)%P;
  ll t=po(A*2%P,P-2,P);
  work(X,((-B-d)*t%P+P)%P);
  work(X,((-B+d)*t%P+P)%P);
}
int main(){
  scanf("%lld",&X);
  if(X>=P)return puts("-1"),0;
  ll A=P-1,B=X,C=X*X%P;
  g[0][0][0]=1;
  g[0][1][1]=1;
  g[1][0][1]=1;
  g[1][1][0]=1;
  g[1][1][1]=1;
  for(int i=2;i<=K;i++)mul(g[i-1],g[1],g[i]);
  w[0][1]=1;
  for(int i=1;i<=K;i++){
    w[i][0]=(g[K][0][0]*w[i-1][0]+g[K][0][1]*w[i-1][1])%P;
    w[i][1]=(g[K][1][0]*w[i-1][0]+g[K][1][1]*w[i-1][1])%P;
  }
  solve(A,B,(C-1+P)%P);
  solve(A,B,(C+1)%P);
  return printf("%lld",ans),0;
}