【BZOJ5104】Fib数列（BSGS，二次剩余）

【BZOJ5104】Fib数列（BSGS，二次剩余）

题面

BZOJ

题解

首先求出斐波那契数列的通项：
令\(A=\frac{1+\sqrt 5}{2},B=\frac{1-\sqrt 5}{2}\)，那么\(f[n]=\frac{1}{\sqrt 5}(A^n-B^n)\)。
然后有\(A=-\frac{1}{B}\)，所以有：\(f[n]=\frac{1}{\sqrt 5}((-\frac{1}{B})^n-B^n)\)。
令\(x=B^n\)，这里需要考虑一下\(n\)的奇偶性：
如果\(n\)是偶数，有：\(\sqrt 5 f[n]=\frac{1}{x}-x\)，转成一元二次方程可以直接求解\(x\)，然后再\(BSGS\)还原。
\(n\)是奇数类似。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
#include<map>
using namespace std;
#define MOD 1000000009
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
namespace SQRT
{
	bool check(int n){return fpow(n,(MOD-1)>>1)==1;}
	int w;struct Num{int a,b;};
	Num operator*(Num a,Num b){return (Num){(1ll*a.a*b.a+1ll*a.b*b.b%MOD*w)%MOD,(1ll*a.a*b.b+1ll*a.b*b.a)%MOD};}
	Num fpow(Num a,int b){Num s=(Num){1,0};while(b){if(b&1)s=s*a;a=a*a;b>>=1;}return s;}
	int Sqrt(int n)
	{
		if(!check(n))return -1;
		int a;do a=rand()%MOD;while(check((1ll*a*a+MOD-n)%MOD));
		w=(1ll*a*a-n+MOD)%MOD;
		return fpow((Num){a,1},(MOD+1)>>1).a;
	}
}
namespace BSGS
{
	map<int,int> M;
	const int m=sqrt(MOD);
	int BSGS(int x,int y)
	{
		M.clear();
		for(int i=0,t=1;i<m;++i,t=1ll*t*x%MOD)M[1ll*y*t%MOD]=i;
		for(int i=1,p=fpow(x,m),t=p;i<=m;++i,t=1ll*t*p%MOD)
			if(M.find(t)!=M.end())
				return i*m-M[t];
		return 2e9;
	}
}
int n,ans=2e9;
int main()
{
	int sqrt5=SQRT::Sqrt(5);
	scanf("%d",&n);
	int A=MOD-1,B=(MOD-1ll*sqrt5*n%MOD)%MOD;
	int d1=(1ll*B*B%MOD+MOD-4ll*A%MOD)%MOD;
	int d2=(1ll*B*B%MOD+4ll*A%MOD)%MOD;
	if(SQRT::check(d1))
	{
		int vd=SQRT::Sqrt(d1);
		int x1=1ll*(0ll+MOD-B+vd)%MOD*fpow(A+A,MOD-2)%MOD;
		int x2=1ll*(0ll+MOD-B+MOD-vd)%MOD*fpow(A+A,MOD-2)%MOD;
		int beta=1ll*(1+MOD-sqrt5)*fpow(2,MOD-2)%MOD;
		int n1=BSGS::BSGS(beta,x1);
		int n2=BSGS::BSGS(beta,x2);
		if(!(n1&1))ans=min(ans,n1);
		if(!(n2&1))ans=min(ans,n2);
	}
	if(SQRT::check(d2))
	{
		int vd=SQRT::Sqrt(d2);
		int x1=1ll*(0ll+MOD-B+vd)%MOD*fpow(A+A,MOD-2)%MOD;
		int x2=1ll*(0ll+MOD-B+MOD-vd)%MOD*fpow(A+A,MOD-2)%MOD;
		int beta=1ll*(1+MOD-sqrt5)*fpow(2,MOD-2)%MOD;
		int n1=BSGS::BSGS(beta,x1);
		int n2=BSGS::BSGS(beta,x2);
		if(n1&1)ans=min(ans,n1);
		if(n2&1)ans=min(ans,n2);
	}
	printf("%d\n",ans<MOD?ans:-1);
	return 0;
}