Fibonacci Number CodeChef - FN (BSGS+cipolla算法)

CodeChef - FN Fibonacci Number(BSGS+cipolla)

CodeChef - FN

题目大意

给出一个数问其最小与斐波那契数列数列的第几项同余于p

解题思路

斐波那契数列有着通项
$$F_n=\frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2}{)}^n-(\frac{1-\sqrt{5}}{2}{)}^n]$$
由此，问题就变成了求解满足
$$F_n\equiv c(modp)$$
的最小n

不妨设$x=\sqrt{5},y=\frac{1+\sqrt{5}}{2}$由此原式就可以转化成
$$y^n-(-y{)}^{-n}\equiv cx(modp)$$
令$q=y^n$并令等式两边同乘上$q$则有
$$q^2-cxq-(-1{)}^n\equiv 0(modp)$$
由求根公式知
$$q=\frac{cx\pm \sqrt{(cx{)}^2+4(-1{)}^n}}{2}$$
由此可以通过cipolla算法先求出q再由BSGS求出y

AC代码

#include<bits/stdc++.h>
using namespace std;
#define int long long 
typedef long long LL;
unordered_map<int,int> mp;
LL quick_pow(LL a,LL b,LL p)
{
	LL ans=1;a%=p;
	while(b)
	{
		if(b&1) ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans;
}
namespace Cipolla{
    LL w,P;
    struct E{
        LL x,y;
        E(LL _x,LL _y):x(_x),y(_y){}
        friend E operator *(E a,E b){
            return E((a.x*b.x%P+a.y*b.y%P*w)%P,(a.x*b.y%P+a.y*b.x%P)%P);
        }
    };
    inline E Pow(E x,int y){
        E ret=E(1,0);
        for(;y;y>>=1,x=x*x) if(y&1) ret=ret*x;
        return ret;
    }
    LL calc(LL x,LL p){
        P=p; x%=P;
        if(x==0) return 0;
        LL tmp=quick_pow(x,(p-1)/2,p);
        if(tmp!=1) return -1;
        while(1){
            tmp=rand()%p;
            LL cur=(tmp*tmp%P+P-x)%P;
            if(quick_pow(cur,(p-1)/2,p)==p-1) break;
        }
        w=(tmp*tmp%P+P-x)%P;
        E A=E(tmp,1);
        return Pow(A,(p+1)/2).x;
    }
}


int BSGS(int B,int N,int P,int sign)//B^L=N(mod p)
{
	mp.clear();
	if(B%P==0) return -1;
	LL ans=0;
	LL m=ceil(sqrt(P));
	for(int i=0;i<=m;i++)
	{
		if(!i) ans=N%P,mp[ans]=i;
		else
		{
			ans=ans*B%P;
			mp[ans]=i;
		}
	}
	LL t=quick_pow(B,m,P);ans=1;
	for(int i=1;i<=m;i++)
	{
		ans=(ans*t)%P;
		if(mp[ans])
		{
			int s=i*m-mp[ans];
			if((s&1)==sign)return (s%P+P)%P;
		}
	}
	return -1;
}
int solve(int c,int p)
{
	int x=Cipolla::calc(5,p);
	int cx=c*x%p;
	int frac12=(p+1)/2;
	int y=(1+x)*frac12%p;
	int ans=p+1;
	int delta=Cipolla::calc(5LL*c*c%p+4,p);
	int yn1=((cx+delta)%p)*frac12%p;
	int yn2=((cx+p-delta)%p)*frac12%p;
	int ason;
	if(~delta){
	ason=BSGS(y,yn1,p,0);
	if(ason!=-1) ans=min(ans,ason);
	ason=BSGS(y,yn2,p,0);
	if(ason!=-1) ans=min(ans,ason);
	}
	delta=Cipolla::calc(5LL*c*c-4,p);
	yn1=((cx+delta)%p)*frac12%p;
	yn2=((cx+p-delta)%p)*frac12%p;
	if(~delta){
	ason=BSGS(y,yn1,p,1);
	if(ason!=-1) ans=min(ans,ason);
	ason=BSGS(y,yn2,p,1);
	if(ason!=-1) ans=min(ans,ason);
	}
	if(ans!=p+1) return ans;
	return -1;
}
int32_t main()
{
	int t;
	cin>>t;
	while(t--)
	{
		int c,p;
		cin>>c>>p;
		cout<<solve(c,p)<<endl;
	}
}