斐波那契數列的通項：
\[\frac{1}{\sqrt{5}}((\frac{1+\sqrt{5}}{2})-(\frac{1-\sqrt{5}}{2}))\]
設T=\(\sqrt{5}*N\),\(y=\frac{\sqrt{5}+1}{2}\)
原式可化為\(y^n-(-\frac{1}{y}^n) \equiv T(mod\ p)\)
我們設\(t=y^n\)
原式可繼續化為\(t-T*t \equiv (-1)^n(mod\ p)\)
然後我們對n進行奇偶討論。
即分別求出\(t-T*t\equiv 1(mod\ p)\)和\(t-T*t\equiv -1(mod\ p)\)的t的解，這個用求根公式+二次剩余求出。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstdlib>
#include<ctime>
using namespace std;
#define int long long
const int p=1e9+9;
const int N=201000;
int head[N],cnt;
int num,ans[10],w;
struct edge{
    int w,id,nxt;
}e[N];
void add(int u,int w,int id){
    for(int i=head[u];i;i=e[i].nxt)
        if(e[i].w==w&&e[i].id<id)swap(e[i].id,id);
    cnt++;
    e[cnt].nxt=head[u];
    e[cnt].w=w;
    e[cnt].id=id;
    head[u]=cnt;
}
struct comple{
    int x,y;
    comple (int xx=0,int yy=0){
        x=xx;y=yy;
    }
};
comple operator *(comple a,comple b){
    return comple(((a.x*b.x%p+a.y*b.y%p*w%p)%p+p)%p,((a.x*b.y%p+a.y*b.x%p)%p+p)%p);
}
int random(int x){
    return rand()*rand()%p;
}
int ksm(int x,int b){
    int tmp=1;
    while(b){
        if(b&1)tmp=tmp*x%p;
        x=x*x%p;
        b>>=1;
    }
    return tmp;
}
comple ksm(comple x,int b){
    comple tmp(1,0);
    while(b){
        if(b&1)tmp=tmp*x;
        x=x*x;
        b>>=1;
    }
    return tmp;
}
int Sqrt(int x){
    if(p==2)return x;
    if(ksm(x,(p-1)/2)+1==p)return -1;
    int a;
    while(233){
        a=random(p);
        w=((a*a%p-x)%p+p)%p;
        if(ksm(w,(p-1)/2)+1==p)break;
    }
    comple res(a,1);
    comple ans(0,0);
    ans=ksm(res,(p+1)/2);
    return ans.x;
}
int BSGS(int a,int b){
    int block=sqrt(p)+1;
    int tmp=b;
    for(int i=0;i<block;i++,tmp=tmp*a%p)add(tmp%200000+1,tmp,i);
    a=ksm(a,block);
    if(a==0)return b==0?1:-1;
    tmp=1;
    for(int i=0;i<=block;i++,tmp=tmp*a%p){
        for(int j=head[tmp%200000+1];j;j=e[j].nxt)
            if(e[j].w==tmp&&i*block-e[j].id>=0)return i*block-e[j].id;
    }
    return -1;
}
int read(){
    int sum=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9'){sum=sum*10+ch-'0';ch=getchar();}
    return sum*f;
}
signed main(){
    srand(time(NULL));
    int n=read()%p;
    int a=Sqrt(5);
    int T=a*n%p;
    int y=(1+a)*ksm(2,p-2)%p;
    int x1=Sqrt(T*T%p+4ll);
    if(x1!=-1){
        int t1=(T+x1)%p*ksm(2,p-2)%p;
        int t2=((T-x1)%p+p)%p*ksm(2,p-2)%p;
        int ans1=BSGS(y,t1);
        cnt=0;memset(head,0,sizeof(head));
        int ans2=BSGS(y,t2);
        if(ans1!=-1)ans[++num]=ans1;
        if(ans2!=-1)ans[++num]=ans2;
    }
    int x2=Sqrt(T*T%p-4);
    if(x2!=-1){
        int t1=(T+x2)%p*ksm(2,p-2)%p;
        int t2=((T-x2)%p+p)%p*ksm(2,p-2)%p;
        cnt=0;memset(head,0,sizeof(head));
        int ans1=BSGS(y,t1);
        cnt=0;memset(head,0,sizeof(head));
        int ans2=BSGS(y,t2);
        if(ans1!=-1)ans[++num]=ans1;
        if(ans2!=-1)ans[++num]=ans2;
    }
    if(num==0)printf("-1");
    else {
        sort(ans+1,ans+1+num);
        printf("%lld",ans[1]);
    }
    return 0;
}
 

BZOJ 5104 Fib數列（二次剩余+BSGS）