【模板】二次剩余Cipolla算法/欧拉准则-bzoj5104: Fib数列

---

欧拉准则

对于质数$p$，$x^2\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}p)\iff a^{\frac{p-1}{2}}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$。

证明：

充分性：
$a^{\frac{p-1}{2}}=(x^2{)}^{\frac{p-1}{2}}=x^{p-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$

必要性：
设$g$为模$p$意义下的原根，$g^i\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}p)$，则
$g^{\frac{i(p-1)}{2}}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}p)$
$g$为原根，则$(p-1)|\frac{i(p-1)}{2}$，$i$为偶数，$x\equiv g^{\frac{i}{2}}(\mathrm{m}\mathrm{o}\mathrm{d}p)$

---

Cipolla算法

一种求解模奇质数意义下的二次同余方程的算法。

即$p$为奇质数，给定$a$，求$x$满足：$$x^2\equiv a(\mathrm{m}\mathrm{o}\mathrm{d}p)$$

复杂度$O(\log p)$

具体算法:

可能需要死记硬背一下。。。

随机找到任意一个$b$，满足$b^2-a\equiv w$，$w$是模$p$意义下的非二次剩余($w^{\frac{p-1}{2}}\equiv -1(\mathrm{m}\mathrm{o}\mathrm{d}p)$)。

类似于虚数，设$i=\sqrt{w}$，扩域后每个数表达为$(a,b)=a+bi$，运算为：
$(r_1,d_1)+(r_2,d_2)=(r_1+r_2,d_1+d_2)$
$(r_1,d_1)\cdot (r_2,d_2)=(r_1{r}_2+d_1{d}_{2}w,r_1{d}_2+r_2{d}_1)$

$<G,+,\cdot >$是个环，满足交换律和结合律。

答案即为：$$x=(b+i{)}^{\frac{p+1}{2}}$$

证明：
$x^2$
$=(b+i{)}^p(b+i)$
$=(b^p+i^p)(b+i)$ ($(b+i{)}^p$二项式定理展开后其余项均含$p$模后为0)
由$b^p\equiv b,i^p=w^{\frac{p-1}{2}}i=-i$
则$x^2=(b-i)(b+i)=b^2-i^2=b^2-w=a$

---

例题

传送门：bzoj5104

求在$(\mathrm{m}\mathrm{o}\mathrm{d}10^9+9)$的意义下，数字$N$在$Fib$数列中出现在哪个位置
无解输出$-1$
$N<=10^9+9$

---

题解

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

吐槽一下这怎么记得到啊。。。

设$t=(\frac{1+\sqrt{5}}{2}{)}^n,T=\sqrt{5}N$

则$t-(-1{)}^n\frac{1}{t}=T$

左右同乘$T$：

得到方程$t^2-Tt-(-1{)}^n=0$

分$n$的奇偶讨论，先求出$\sqrt{\Delta }$后$BSGS$求解。

---

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+9,qrt5=383008016,nv2=500000005,blk=sqrt(mod)+1,bs=691504013;
int n,ans,ca,cb;
map<int,int>mp[2];

inline int fp(int x,int y)
{
	int re=1;
	for(;y;y>>=1,x=(ll)x*x%mod)
	  if(y&1) re=(ll)re*x%mod;
	return re;
}

inline int ad(int x,int y){x+=y;return x>=mod?x-mod:x;}
inline int dc(int x,int y){x-=y;return x<0?x+mod:x;}

namespace sc_surplus{
    int w,a,b;
    
    struct cc{
    	int r,i;
    	cc(int r,int i):r(r),i(i){};
    	cc operator *(const cc&ky)const
		{return cc(ad((ll)r*ky.r%mod,(ll)i*ky.i%mod*(ll)w%mod),ad((ll)r*ky.i%mod,(ll)i*ky.r%mod));}
        cc operator *=(const cc&ky){return (*this)=(*this)*ky;}
	};
    
    inline int fp(cc x,int y)
    {
    	cc re=cc(1,0);
    	for(;y;y>>=1,x*=x) if(y&1) re*=x;
    	return re.r;
    }
    
    inline bool ck(int x)
    {return (::fp((w=dc((ll)x*x%mod,a)),(mod-1)/2)==mod-1);}
    
	inline int fd_sqrt(int x)
	{
		if(::fp(x,(mod-1)/2)!=1) return 0;
		a=x;for(b=rand();!ck(b);b=rand());
		return fp(cc(b,1),(mod+1)/2);
	}	
}

inline void cal(int od)
{
	int i,j,a=ca,b=cb,pw;
	for(i=0;i<2;++i) if(!mp[i].empty()) mp[i].clear();
	for(i=0;i<blk;++i){
		j=i&1;
		if(mp[j].find(a)==mp[j].end()) mp[j][a]=i;
		if(mp[j].find(b)==mp[j].end()) mp[j][b]=i;
		a=(ll)a*bs%mod;b=(ll)b*bs%mod;
	}
	pw=fp(bs,blk);int vl=pw,re=blk;
	for(i=0;i<=blk;++i){
		j=-1;
		if(od){
		   if((re&1)&&(mp[0].find(vl)!=mp[0].end())) j=mp[0][vl];
		   else if((!(re&1))&&(mp[1].find(vl)!=mp[1].end())) j=mp[1][vl];
		}else{
		   if((re&1)&&(mp[1].find(vl)!=mp[1].end())) j=mp[1][vl];
		   else if((!(re&1))&&(mp[0].find(vl)!=mp[0].end())) j=mp[0][vl];
		}if(~j) {re-=j;ans=min(ans,re);return;}
		vl=(ll)vl*pw%mod;re+=blk;
	}
}

int main(){
	srand(time(0));
	scanf("%d",&n);ans=mod+10;
	if(n==1) {puts("1");return 0;}
	int i,j,x,y;
	n=(ll)n*qrt5%mod;
	if((x=sc_surplus::fd_sqrt(((ll)n*n+4)%mod))){
		ca=(ll)ad(n,x)*nv2%mod;cb=(ll)dc(n,x)*nv2%mod;cal(0);
	}
	if((x=sc_surplus::fd_sqrt(((ll)n*n-4+mod)%mod))){
		ca=(ll)ad(n,x)*nv2%mod;cb=(ll)dc(n,x)*nv2%mod;cal(1);
	}
	printf("%d",ans==mod+10?(-1):ans);
	return 0;
}