【模闆】二次剩餘Cipolla算法/歐拉準則-bzoj5104: Fib數列

歐拉準則

對于質數 p p p， x 2 ≡ a ( m o d p ) ⇔ a p − 1 2 ≡ 1 ( m o d p ) x^2\equiv a\pmod p\Leftrightarrow a^{\frac{p-1}{2}}\equiv 1\pmod p x2≡a(modp)⇔a2p−1​≡1(modp)。

證明：

充分性：

a p − 1 2 = ( x 2 ) p − 1 2 = x p − 1 ≡ 1 ( m o d p ) a^{\frac{p-1}{2}}=(x^2)^{\frac{p-1}{2}}=x^{p-1}\equiv 1\pmod p a2p−1​=(x2)2p−1​=xp−1≡1(modp)

必要性：

設 g g g為模 p p p意義下的原根， g i ≡ a ( m o d p ) g^i\equiv a\pmod p gi≡a(modp)，則

g i ( p − 1 ) 2 ≡ 1 ( m o d p ) g^{\frac{i(p-1)}{2}}\equiv 1\pmod p g2i(p−1)​≡1(modp)

g g g為原根，則 ( p − 1 ) ∣ i ( p − 1 ) 2 (p-1)|\frac{i(p-1)}{2} (p−1)∣2i(p−1)​， i i i為偶數， x ≡ g i 2 ( m o d p ) x\equiv g^{\frac i2}\pmod p x≡g2i​(modp)

Cipolla算法

一種求解模奇質數意義下的二次同餘方程的算法。

即 p p p為奇質數，給定 a a a，求 x x x滿足： x 2 ≡ a ( m o d p ) x^2\equiv a\pmod p x2≡a(modp)

複雜度 O ( log ⁡ p ) O(\log p) O(logp)

具體算法:

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

随機找到任意一個 b b b，滿足 b 2 − a ≡ w b^2-a\equiv w b2−a≡w， w w w是模 p p p意義下的非二次剩餘( w p − 1 2 ≡ − 1 ( m o d p ) w^{\frac{p-1}{2} }\equiv -1\pmod p w2p−1​≡−1(modp))。

類似于虛數，設 i = w i= \sqrt w i=w

​，擴域後每個數表達為 ( a , b ) = a + b i (a,b)=a+bi (a,b)=a+bi，運算為：

( r 1 , d 1 ) + ( r 2 , d 2 ) = ( r 1 + r 2 , d 1 + d 2 ) (r_1,d_1)+(r_2,d_2)=(r_1+r_2,d_1+d_2) (r1​,d1​)+(r2​,d2​)=(r1​+r2​,d1​+d2​)

( r 1 , d 1 ) ⋅ ( r 2 , d 2 ) = ( r 1 r 2 + d 1 d 2 w , r 1 d 2 + r 2 d 1 ) (r_1,d_1)·(r_2,d_2)=(r_1 r_2+d_1d_2w,r_1d_2+r_2d_1) (r1​,d1​)⋅(r2​,d2​)=(r1​r2​+d1​d2​w,r1​d2​+r2​d1​)

&lt; G , + , ⋅ &gt; &lt;G,+,·&gt; <G,+,⋅>是個環，滿足交換律和結合律。

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

證明：

x 2 \quad x^2 x2

= ( b + i ) p ( b + i ) =(b+i)^p(b+i) =(b+i)p(b+i)

= ( b p + i p ) ( b + i ) =(b^p+i^p)(b+i) =(bp+ip)(b+i) ( ( b + i ) p (b+i)^p (b+i)p二項式定理展開後其餘項均含 p p p模後為0)

由 b p ≡ b , i p = w p − 1 2 i = − i b^p\equiv b,i^p=w^{\frac{p-1}{2}}i=-i bp≡b,ip=w2p−1​i=−i

則 x 2 = ( b − i ) ( b + i ) = b 2 − i 2 = b 2 − w = a x^2=(b-i)(b+i)=b^2-i^2=b^2-w=a x2=(b−i)(b+i)=b2−i2=b2−w=a

例題

傳送門：bzoj5104

求在 ( m o d 1 0 9 + 9 ) \pmod{10^9+9} (mod109+9)的意義下，數字 N N N在 F i b Fib Fib數列中出現在哪個位置

無解輸出 − 1 -1 −1

N &lt; = 1 0 9 + 9 N &lt; = 10^9+9 N<=109+9

題解

F i b Fib Fib數列的通項公式： F n = ( 1 + 5 2 ) n − ( 1 − 5 2 ) n 5 F_n=\frac{(\frac{1+\sqrt 5}{2})^n-(\frac{1-\sqrt 5}{2})^n}{\sqrt 5} Fn​=5

​(21+5

​​)n−(21−5

​​)n​

吐槽一下這怎麼記得到啊。。。

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

​​)n,T=5

​N

則 t − ( − 1 ) n 1 t = T t-(-1)^n\frac{1}{t}=T t−(−1)nt1​=T

左右同乘 T T T：

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

分 n n n的奇偶讨論，先求出 Δ \sqrt\Delta Δ

​後 B S G S BSGS 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;
}