BZOJ5104：Fib数列（二次剩余+BSGS）

传送门

题意：
给定$a$，求$F_n \equiv a \pmod{1e9+9}$，输出最小的$n$。

题解：
Fib数列的通项：$F_n=\frac{1}{\sqrt{5}}((\frac{1+\sqrt{5}}{2})^x-\frac{1-\sqrt{5}}{2})^n)$。

$5$是$1e9+9$的一类二次剩余。可以对其开根号：

$$5\equiv 383008016^2\pmod{1e9+9}$$
（开根号的方法： http://blog.csdn.net/qq_35649707/article/details/78922508）

设$x=\frac{1+\sqrt{5}}{2}$，$T=\sqrt{5}a$那么公式可表达为:

$$x^n-(-\frac{1}{x})^n \equiv T\pmod{1e9+9}$$ 。

设$t=x^n$，则：

$$t-(-1)^n \frac{1}{t}\equiv T \pmod{1e9+9}$$

不妨奇偶讨论，两边同时乘上$t$，可得：

$$\begin{cases} t^2+ 1\equiv T*t (n \equiv 1 \pmod 2 )\\ t^2-1\equiv T*t (n \equiv 0 \pmod 2) \end{cases}$$

如果得到了这个方程的$t$的解，那么可以套用$BSGS$解出答案，至于解这个方程，只需套用求根公式。解$\sqrt{\Delta }$再次套用二次剩余即可。

#include<bits/stdc++.h>
using namespace std;
const int Mod=1e9+9;
const int Sqrt5=383008016,InvS=276601605,Inv2=500000005,A=691504013;
inline int power(int a,int b){
    int rs=1;
    for(;b;b>>=1,a=1ll*a*a%Mod)if(b&1)rs=1ll*rs*a%Mod;
    return rs;
}
int n,S,t,t1,t2;
namespace findroot{
    int w=0,a,b;
    struct Complex{
        int r,i;
        Complex(int r,int i):r(r),i(i){}
        friend inline Complex operator *(const Complex &a,const Complex &b){
            return Complex((1ll*a.r*b.r+1ll*a.i*b.i%Mod*w%Mod)%Mod,(1ll*a.r*b.i+1ll*a.i*b.r)%Mod);
        }
    };
    inline bool check(int v){
        w=(1ll*v*v-a+Mod)%Mod;
        return ::power(w,(Mod-1)/2)!=1;
    }
    inline int power(Complex a,int b){
        Complex rs(1,0);
        for(;b;b>>=1,a=a*a)if(b&1)rs=rs*a;
        return rs.r;
    }
    inline int calc(int k){
        a=k;
        if(::power(a,(Mod-1)/2)!=1)return 0;
        while(!check(b=rand()));
        return power(Complex(b,1),(Mod+1)/2);
    }
}
const int H=400009,Bl=sqrt(Mod)+1;
struct Hash_Table{
    vector< pair<int,int> >vv[H];
    inline void insert(int p,int i){
        vv[p%H].push_back(make_pair(p,i));
    }
    inline int query(int p){
        int v=p%H,mx=-1;
        for(int i=vv[v].size()-1;i>=0;i--)
            if(vv[v][i].first==p)mx=max(mx,vv[v][i].second);
        return mx;
    }
}h1,h2,h3,h4;
inline long long findans(int bz,int rt,Hash_Table &hh1,Hash_Table&hh2){
    int pw1=1;
    for(int i=0;i<Bl;i++){
        if(i&1){
            hh1.insert(1ll*pw1*t1%Mod,i);
            hh1.insert(1ll*pw1*t2%Mod,i);
        }else{
            hh2.insert(1ll*pw1*t1%Mod,i);
            hh2.insert(1ll*pw1*t2%Mod,i);
        }
        pw1=1ll*pw1*rt%Mod;
    }
    int now1=pw1;long long pos=Bl;
    do{
        int p1=(bz)?((pos&1)?hh2.query(now1):hh1.query(now1)):((pos&1)?hh1.query(now1):hh2.query(now1));
        if(p1!=-1){
            pos=pos-p1;
            break;
        }
        now1=1ll*now1*pw1%Mod;
        pos+=Bl;
    }while(pos<=INT_MAX);
    return pos;
}
int main(){
    srand(23333);
    scanf("%d",&n);
    if(n==1){puts("1");return 0;}
    n=1ll*n*Sqrt5%Mod;
    int rt1=findroot::calc((1ll*n*n+4)%Mod);
    int rt2=findroot::calc((1ll*n*n-4)%Mod);
    long long ans=INT_MAX;
    if(rt1){
        t1=1ll*(n+rt1)%Mod*Inv2%Mod;
        t2=1ll*(n-rt1+Mod)%Mod*Inv2%Mod;
        ans=min((long long)ans,findans(0,A,h1,h2));
    }
    if(rt2){
        t1=1ll*(n+rt2)%Mod*Inv2%Mod;
        t2=1ll*(n-rt2+Mod)%Mod*Inv2%Mod;
        ans=min((long long)ans,findans(1,A,h3,h4));
    }
    printf("%d\n",(ans<INT_MAX?ans:-1));
}