數論之BSGS算法

UPDATE ON 2020.12.23 更新了原代碼中的錯誤

BSGS算法

\(BSGS(baby-step giant-step)\)，即大步小步算法。經常使用於求解離散對數問題。

形式化地說，該算法能夠在 \(O \sqrt{p}\) 的時間內求解

\[a^x \equiv b (mod\ p) \]

其中 \(gcd(a,p)=1\)

根據歐拉定理，當 \(gcd(a,p)=1\) 時，\(a^{\varphi(p)}mod\ p=1\)

而 \(a^0mod\ p=1\)

因此從 \(1 \sim \varphi(p)\) 出現了一個循環節

所以，咱們只須要求出 \(a^1 \sim a^{\varphi(p)}\)

而後分別帶入左右兩邊判斷是否相等

通常的題目中 \(p\) 爲質數，因此 \(\varphi(p)=p-1\)

通常就用 \(p\) 替代

可是這樣作時間複雜度過高，能夠用分塊的思想優化

令 \(a^x=a^{km-t}\)，那麼原式就變成了

\[a^{km} \equiv ba^t (mod\ p) \]

咱們只須要分別求出 \(t=1 \sim m\) 時 \(ba^t\) 的值，將其存到哈希表裏

而後從 \(1 \sim p/m\) 枚舉 \(k\) ，判斷哈希表中有沒有和 \(a^{km}\) 相等的數

當 \(m=\sqrt{p}\) 時最優

代碼

模板題

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#define rg register
inline int read(){
   rg int x=0,fh=1;
   rg char ch=getchar();
   while(ch<'0' || ch>'9'){
   	if(ch=='-') fh=-1;
   	ch=getchar();
   }
   while(ch>='0' && ch<='9'){
   	x=(x<<1)+(x<<3)+(ch^48);
   	ch=getchar();
   }
   return x*fh;
}
const int maxn=1e5+5,mod=1e5+3;
int p,a,n,blo,h[maxn],tot=1;
struct asd{
   int nxt,val,num;
}b[maxn];
void ad(int val,int id){
   rg int now=val%mod;
   for(rg int i=h[now];i!=-1;i=b[i].nxt){
   	rg int cs=b[i].val;
   	if(cs==val){
   		if(b[i].num<id) b[i].num=id;
   		return;
   	}
   }
   b[tot].val=val;
   b[tot].nxt=h[now];
   b[tot].num=id;
   h[now]=tot++;
}
int cx(int val){
   rg int now=val%mod;
   for(rg int i=h[now];i!=-1;i=b[i].nxt){
   	rg int cs=b[i].val;
   	if(cs==val) return b[i].num;
   }
   return -1;
}
int main(){
   memset(h,-1,sizeof(h));
   p=read(),a=read(),n=read();
   blo=sqrt(p)+1;//不要忘了加上一個1保證正確性
   rg int now=n,cs=1,haha;
   for(rg int i=1;i<=blo;i++){
   	now=1LL*now*a%p;
   	ad(now,i);
   	cs=1LL*cs*a%p;
   }
   now=1;
   for(rg int i=1;i<=blo;i++){
   	now=1LL*now*cs%p;
   	haha=cx(now);
   	if(haha!=-1){
   		printf("%d\n",i*blo-haha);
   		return 0;
   	}
   }
   printf("no solution\n");
   return 0;
}

擴展BSGS算法

求解

\[a^x \equiv b (mod\ p) \]

其中 \(a,p\) 不必定互質

這個時候就不能再找 \(1 \sim \varphi(p)\) 的循環節了，由於不知足歐拉定理的條件

咱們把式子變一下，能夠寫成

\[a \times a^{x-1}+kp=b \]

的形式

根據裴蜀定理，當 \(gcd(a,p)|b\) 時方程有解，不然無解

因此在有解的狀況下，咱們能夠把方程左右兩邊同時除以 \(gcd(a,p)\)，等號不變

此時原方程變爲

\[\frac{a}{gcd(a,p)} \times a^{x-1}+k \frac{p}{gcd(a,p)}=\frac{b}{gcd(a,p)} \]

把 \(\frac{p}{gcd(a,p)}\) 當作新的 \(p\)

把 \(\frac{b}{gcd(a,p)}\) 當作新的 \(b\)

同時記錄一下左邊的係數之積 \(\frac{a}{gcd(a,p)}\)，並不斷地提出 \(a\)

這樣在最後必定會出現 \(gcd(a,p)=1\) 的狀況，用普通的 \(bsgs\) 便可

不要忘了把提出去的 \(a\) 的係數加上

代碼

模板題

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#define rg register
inline int read(){
	rg int x=0,fh=1;
	rg char ch=getchar();
	while(ch<'0' || ch>'9'){
		if(ch=='-') fh=-1;
		ch=getchar();
	}
	while(ch>='0' && ch<='9'){
		x=(x<<1)+(x<<3)+(ch^48);
		ch=getchar();
	}
	return x*fh;
}
const int maxn=1e5+5,mod=1e5+3;
int h[maxn],tot=1;
struct asd{
	int nxt,val,num;
}b[maxn];
void ad(int val,int id){
	rg int now=val%mod;
	for(rg int i=h[now];i!=-1;i=b[i].nxt){
		rg int cs=b[i].val;
		if(cs==val){
			if(b[i].num<id) b[i].num=id;
			return;
		}
	}
	b[tot].val=val;
	b[tot].nxt=h[now];
	b[tot].num=id;
	h[now]=tot++;
}
int cx(int val){
	rg int now=val%mod;
	for(rg int i=h[now];i!=-1;i=b[i].nxt){
		rg int cs=b[i].val;
		if(cs==val) return b[i].num;
	}
	return -1;
}
int bsgs(rg int a,rg int p,rg int n,rg int xs){
	memset(h,-1,sizeof(h));
	tot=1;
	rg int blo=sqrt(p)+1;
	rg int ac1=n,ac2=1,ac3=xs;
	for(rg int i=1;i<=blo;i++){
		ac1=1LL*ac1*a%p;
		ad(ac1,i);
		ac2=1LL*ac2*a%p;
	}
	for(rg int i=1;i<=blo;i++){
		ac3=1LL*ac3*ac2%p;
		ac1=cx(ac3);
		if(ac1!=-1) return i*blo-ac1;
	}
	return -1;
}
int gcd(rg int aa,rg int bb){
	return bb==0?aa:gcd(bb,aa%bb);
}
int main(){
	rg int a,p,n,haha=1,now,ans,cnt;
	rg bool jud=0;
	while(1){
		a=read(),p=read(),n=read();
		if(a==0 && p==0 && n==0) break;
		haha=1,jud=0,cnt=0;
		a%=p,n%=p;
		if(n==1){
			printf("0\n");
		} else if(a==0 && n==0){
			printf("0\n");
		} else if(a==0){
			printf("No Solution\n");
		} else if(n==0){
			while(1){
				if(p==1){
					printf("%d\n",cnt);
					break;
				}
				now=gcd(a,p);
				if(now==1 && p!=-1){
					printf("No Solution\n");
					break;
				}
				cnt++,p/=now;
			}
		} else {
			while(1){
				now=gcd(a,p);
				if(n%now!=0){
					jud=1;
					break;
				}
				if(now==1){
					ans=bsgs(a,p,n,haha);
					if(ans==-1) jud=1;
					ans+=cnt;
					break;
				}
				cnt++,p/=now,n/=now;
				haha=1LL*haha*(a/now)%p;
				if(haha==n){
					ans=cnt;
					break;
				}
			}
			if(jud){
				printf("No Solution\n");
			} else {
				printf("%d\n",ans);
			}
		}
	}
	return 0;
}