【原创】BSGS和EXBSGS

气拔山兮力盖世

说在前面

数论不知道是复习还是学习了，Part 5。

BSGS

BSGS算法就是用来已知$y,z,p$，求解$y^x\equiv z\mathrm{m}\mathrm{o}\mathrm{d}p$的最小自然数解的，暂时要求$y,p$互质。

取一个数$m$，然后我们用数$a,b$使得$x=am-b$，自然$b\in [0,m-1]$。
则原方程$⟺y^{am}\equiv z{y}^b\mathrm{m}\mathrm{o}\mathrm{d}p$。
我们枚举$b\in [0,m-1]$，求出这$m-1个z{y}^b\mathrm{m}\mathrm{o}\mathrm{d}p$的值，存起来。
然后我们枚举$a$，看存起来的桶子里面$y^{ax}$有没有出现过，出现过了我们就找到了一个$x=am-b$了。
如果$am>p$了，就停下来，就没找着。（由费马小定理，再往下就会循环）

然后时间复杂度易知$m=\sqrt{p}$的时候最佳，为$\Theta (\sqrt{p})$。

	if(gcd(y,p)!=1) puts("Orz, I cannot find x!");
	else if(z%p==0) puts("0");
	else
	{
		bool ok=0;
		y%=p,z%=p,M.clear();
		int bas=sqrt(p)+1,tim=ksm(y,bas,p);
		for(int i=0,now=z;i<=bas-1;i++,now=1ll*now*y%p) M[now]=i;
		for(int i=0,now=1;i<=bas+1 && !ok;i++,now=1ll*now*tim%p)
			if(M.count(now) && i*bas-M[now]>0) printf("%d\n",i*bas-M[now]),ok=1;
		if(!ok) puts("Orz, I cannot find x!");
	}

代码是SDOI 2242计算器里抠出来的，用的map所以效率较低，没开O2总共跑了1.89s。
可以邻接表挂链。
（但我没写过）

拓展BSGS

如果$y,p$不互质，但$y$不能是$p$的倍数（废话）。
设$d_1=gcd(y,p)$，如果$d_1\nmid z$，则无解。
否则，两边同时除以$d_1$。
$y^{x-1}\ast \frac{y}{d_1}\equiv \frac{z}{d_1}\mathrm{m}\mathrm{o}\mathrm{d}\frac{p}{d_1}$

设$d_2=gcd(y,\frac{p}{d_1})$，如果$d_2\nmid \frac{z}{d_1}$，则无解。
如果$d_2=1$等会儿再说，否则两边同时除以$d_2$，得到：$y^{x-2}\ast \frac{y^2}{d_1\ast d_2}\equiv \frac{z}{d_1\ast d_2}\mathrm{m}\mathrm{o}\mathrm{d}\frac{p}{d_1\ast d_2}$

重复这个操作k次，直到$d_{k+1}=1$。
此时有$y^{x-k}\ast \frac{y^k}{d_1\ast d_2\ast \cdots \ast d_k}\equiv \frac{z}{{\prod }_{i=1}^k{d}_i}\mathrm{m}\mathrm{o}\mathrm{d}\frac{p}{\prod k}$。
令$c=\frac{y^k}{\prod k},z^{\prime }=\frac{z}{\prod k},p^{\prime }=\frac{p}{\prod k}$。

则原式转化为$y^{x-k}c\equiv z^{\prime }\mathrm{m}\mathrm{o}\mathrm{d}{p}^{\prime }$。
因为$gcd(c,p^{\prime })=1$（肯定的，都除了那么多次了），所以我们可以找到$c$关于$p^{\prime }$的逆元$a$，在两边同时乘以$a$得$y^{x-k}\equiv z^{\prime }a\mathrm{m}\mathrm{o}\mathrm{d}{p}^{\prime }$。其中$(y,p^{\prime })$互质。

这是什么？这是BSGS模板题，快切掉。
注意我们求出的是$x-k$，最后求得的$x$总是不小于$k$的，于是我们还要枚举$x\in [0,k-1]$。
P.S.这一条你可以在还在算k的时候就做了。

代码是洛谷的P4195 【模板】exBSGS/Spoj3105 Mod的。

刚刚我说不写哈希表挂链的时候，我就心底一沉，就有了一种map会T的感觉。
果然。
开O2或者用unordered_map能过，不过为了自己着想还是照着yyb敲了一发hash。
信竞黄旭东就是ko no CRLOSS哒。

当然也有一些显然的优化啦，exgcd就是可以省略的。
但是算了，管他的。

代码

#include<map>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;

inline void Read(int &p)
{
	p=0;
	char c=getchar();
	while(c<'0' || c>'9') c=getchar();
	while(c>='0' && c<='9')
		p=p*10+c-'0',c=getchar();
}

const int mod=100007;
struct hash_table
{
	int head[mod+5],cnt,pos[mod+5],val[mod+5],nxt[mod+5];
	inline void addedge(int u,int v,int w){nxt[++cnt]=head[u],head[u]=cnt,pos[cnt]=w,val[cnt]=v;}
	inline void clear(){memset(head,0,sizeof head),cnt=0;}
	inline void insert(int x,int k){addedge(x%mod,k,x);}
	int count(int x)
	{
		for(int i=head[x%mod];i;i=nxt[i])
			if(pos[i]==x) return val[i];
		return -1;
	}
}M;
int y,z,p,ans;

int gcd(int a,int b){return b?gcd(b,a%b):a;}

int ksm(int a,int b,int p)
{
	int ans=1; a%=p;
	while(b)
	{
		if(b&1) ans=1ll*ans*a%p;
		a=1ll*a*a%p,b>>=1;
	}
	return ans;
}

int bsgs(int y,int z,int p)
{
	if(gcd(y,p)!=1) return -0x3f3f3f3f;
	if(z%p==1) return 0;
		
	y%=p,z%=p,M.clear();
	int bas=sqrt(p)+1,tim=ksm(y,bas,p);
	for(int i=0,now=z;i<=bas-1;i++,now=1ll*now*y%p) M.insert(now,i);
	for(int i=0,now=1;i<=bas+1;i++,now=1ll*now*tim%p)
	{
		int orz=M.count(now);
		if(orz!=-1 && i*bas-orz>=0) return i*bas-orz;
	}
	return -0x3f3f3f3f; 
}

int exgcd(int a,int b,int &x,int &y)
{
	if(!b) {x=1,y=0; return a;}
	int w=exgcd(b,a%b,y,x);
	y-=1ll*(a/b)*x;
	return w;
}

int exbsgs(int y,int z,int p)
{
	int k=0,d,c=1,x0,y0; //c=(y^k)/(d1*d2*d3*...*dk)
	y%=p,z%=p; 
	if(z==1) return 0;
	while((d=gcd(y,p))!=1)
	{
		if(z%d) return -1;
		k++,z/=d,p/=d,c=1ll*c*(y/d)%p;
		if(c==z) return k;
	}
	if(p==1) return k;
	exgcd(c,p,x0,y0),z=(1ll*z*x0%p+p)%p,y%=p;
	return bsgs(y,z,p)+k;
}

int main()
{
	while(1)
	{
		Read(y),Read(p),Read(z);
		if(!y && !z && !p) break;
		ans=exbsgs(y,z,p);
		if(ans>=0) printf("%d\n",ans);
		else puts("No Solution");
	}
}

参考文献

orz yyb 你看我给的是$y,z,p$就知道我是从$\mathcal{Y}\mathcal{Y}\mathcal{B}$那里来的了
orz Miracle
orz DQY

说在后面

我太弱了