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

斐波那契数列的通项：

\[\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的解，这个用求根公式+二次剩余求出。
最后离散对数求出n。
（我写的时候求根公式背错了调了半个小时。。）

#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;
}