高次同余方程（BSGS算法）

---

一、朴素BSGS

1.适用情况

常用来解决形如$b^x\equiv n(mod$ $p)$的问题，且p为质数，b与p互质。

首先根据费马小定理，此时存在式子：$b^{p-1}\equiv 1(mod$ $p)$
而x可以表示为$k\ast (p-1)+t$的形式，说明了此时存在一个循环节，通过枚举$[0,p-1]$内的数可以得知是否存在解。

当p较小的时候可以通过枚举来实现，但如果p过大，枚举的时间复杂度就太高了，这个时候就需要使用朴素BSGS，可以在$O(\sqrt{p})$的时间内完成

2.具体实现

原问题的形式为$b^x\equiv n(mod$ $p)$

先令$m=ceil(\sqrt{p})$，对其进行向上取整，保证后面的操作不会漏掉边界。
那么如果存在解x，x可以写为$i\ast m+j$ $(m>j\ge 0)$

问题就变为了：$b^{i\ast m+j}\equiv n(mod$ $p)$
求是否存在合理的$i,j$

1. 第一种处理思路：
   通过同余式计算，可以通过消去律先把式子变为：$b^j\equiv n\ast b^{-(i\ast m)}(mod$ $\frac{p}{gcd(b^{i\ast m},p)})$
   但其中存在逆元，会让操作变得复杂。

2. 第二种处理思路：
   将x表示为$I\ast m-j$ $(m>j\ge 0)$,这个$I$与上方的$i$相比较，$I=i+1$,那么式子左右同时乘以$b^j$，可以写为:$b^{I\ast m}\equiv n\ast b^j(mod$ $p)$
   这样的话，我们接下来操作就变得简单了。

在同余式：$b^{I\ast m}\equiv n\ast b^j(mod$ $p)$中
先进行$BabyStep$：枚举$j$，把结果存入map里面

	for (int i=0;i<=m;i++){//baby step 
       		if (i==0){
       	   		ans=n%p; mp[ans]=i; continue;
         	}
        	ans=(ans*b)%p;
        	mp[ans]=i;
    	}

再进行$GiantStep$：从小到大枚举$I$，如果找到一个$b^I\%p$在map里面有对应的$j$，那么我们就找到了最小的一个解，此时$x=I\ast m-j$可以计算出来。

		ll t=Pow(b,m,p);
		ans=1;
    	for (int i=1;i<=m;i++){//gaint step
        	ans=(ans*t)%p;
        	if(mp[ans]){
        		int t=i*m-mp[ans];
            	printf("%d\n",(t%p+p)%p);
            	s=true;
            	break; 
         	}
    	}

如果不存在这么一个$I$的话，说明无解。

---

二、例题·Discrete Logging

Discrete Logging

模板题

代码：

#include<cstdio>    
#include<map>  
#include<cmath>   
using namespace std; 
typedef long long ll; 
ll p,b,n;  
map<ll,int> mp;  
 
ll Pow(ll x,ll y,ll m){
	ll ans=1;
	while(y){
		if(y&1)
			ans=(ans*x)%m;
		y>>=1;
		x=(x*x)%m;//更新底数 
	}
	return ans%m;
}
ll bsgs(ll p,ll b,ll n){//b^x=n(mod p)
	mp.clear();

	ll m=ceil(sqrt(p));
	ll ans;
	for (int i=0;i<=m;i++){//baby step 
       		if (i==0){
       	   		ans=n%p; mp[ans]=i; continue;
         	}
        	ans=(ans*b)%p;
        	mp[ans]=i;
    	}   	
		ll t=Pow(b,m,p);
		ans=1;
    	for (int i=1;i<=m;i++){//gaint step
        	ans=(ans*t)%p;
        	if(mp[ans]){
        		int t=i*m-mp[ans];
            	return(t%p+p)%p;
         	}
    	}
    	return -1;
}
int main(){
	while (scanf("%lld%lld%lld",&p,&b,&n)!=EOF){  // b^x=n(mod p) 
		
		ll tmp=bsgs(p,b,n);
		if(tmp==-1)printf("no solution\n");
		else printf("%lld\n",tmp);
	}
}

---

三、拓展BSGS

1.区别

仍然是解决：$b^x\equiv n(mod$ $p)$，但此时$b,p$不互质

2.具体实现

我们的目标是对原式进行处理使之可以使用之前的BSGS算法，但之前的算法有个限定条件：$b，p$互质。

我们先把式子写为：$b^x=y\ast p+n$
$d_1=gcd(p,b)$
再左右同时除以$d_1$：$b^{x-1}\frac{b}{d_1}=y\ast \frac{p}{d_1}+\frac{n}{d_1}$
若$n\%d_1!=0$，则无解；
若$gcd(p,b)！=1$，则重复该操作
……
可得式子：$b^{x-k}\frac{b^k}{d_1\dots d_k}=y\ast \frac{p}{d_1\dots d_k}+\frac{n}{d_1\dots d_k}$
化为同余式即：$b^{x-k}\frac{b^k}{d_1\dots d_k}\equiv \frac{n}{d_1\dots d_k}(mod$ $\frac{p}{d_1\dots d_k})$

不妨以$b^{{}^{\prime }}$表示$b^{x-k}$，以$K$表示$\frac{b^k}{d_1\dots d_k}$，$n^{{}^{\prime }}$表示$\frac{n}{d_1\dots d_k}$，以$p^{{}^{\prime }}$表示$\frac{p}{d_1\dots d_k})$

得到同余式的另一个表达：$b^{{}^{\prime }}K\equiv n^{{}^{\prime }}(mod$ $p^{{}^{\prime }})$
由于$K$与$p^{{}^{\prime }}$互质，可以直接使用消去律，得：$b^{{}^{\prime }}\equiv K^{-1}\ast n^{{}^{\prime }}(mod$ $p^{{}^{\prime }})$

此时的$b^{{}^{\prime }}，K^{-1}\ast n^{{}^{\prime }}，p^{{}^{\prime }}$为常数，可以直接使用上面的思路，进行处理。

需要注意的是，此时得到的解不一定是最小的解，还需要验证一下在当前结果之前有没有更小的解，不过这个验证一般不会太耗时。

---

四、例题·Clever Y

Clever Y

此时如果使用map，将会超时，所以需要自己使用哈希表来实现功能。

#include<cstdio>    
#include<map>  
#include<cmath>  
#include<cstring> 
using namespace std; 
typedef long long ll; 
ll p,b,n;  
const int HASH = 65536;
 
struct hashMap
{
    int head[HASH], Size;
    int next[HASH];
    ll state[HASH], val[HASH];
 
    void init()
    {
        memset(head, -1, sizeof(head));
        Size = 0;
    }
    void Insert(ll st, ll sv)
    {
        ll h = st % HASH;
        for (int p = head[h]; p != -1; p = next[p])
            if (state[p] == st)
                return;
        state[Size] = st;
        val[Size] = sv;
        next[Size] = head[h];
        head[h] = Size++;
    }
    ll Find(ll st)
    {
        ll h = st % HASH;
        for (int p = head[h]; p != -1; p = next[p])
            if (state[p] == st)
                return val[p];
        return -1;
    }
} hashmap;

ll Pow(ll x,ll y,ll m){
	ll ans=1;
	while(y){
		if(y&1)
			ans=(ans*x)%m;
		y>>=1;
		x=(x*x)%m;//更新底数 
	}
	return ans%m;
}

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

ll exgcd(ll a, ll b, ll &x, ll &y){//最终返回的是gcd(a,b)
	if(!b) {
		x=1;y=0;
		return a; //到达递归边界返回
	}
	ll r=exgcd(b,a%b,x,y);
	ll temp=y; //把x y变成上一层的
	y=x-a/b*y;
	x=temp;
	return r; //得到a b的最大公因数
}

ll Inval(ll a,ll n){//求a模n的逆元 
    ll x,y;
    ll d=exgcd(a,n,x,y);
    return d==1?(x%n+n)%n:-1;
}

ll bsgs(ll p,ll b,ll n){//b^x=n(mod p)
	
	for(ll i=0,e=1;i<=100;i++){//由于数据大小，可以先验证有没有解
		if(e==n){
			return i;
		}
		e=e*b%p;
	}
		
	ll k=1,cnt=0;
	while(1){
		ll tmp=gcd(p,b);
		if(tmp==1)break;
		if(n%tmp!=0){
			return -1;
		}
		n/=tmp;
		p/=tmp;
		k=(k*b/tmp)%p;
		cnt++;
	}
		
	ll m=ceil(sqrt(p));
	ll ans=1;
	hashmap.init();
	
	ll K=Inval(k,p); 
	for (int i=0;i<=m;i++){//baby step 
       		if (i==0){
       	   		ans=K*n%p; hashmap.Insert(ans,i); continue;
         	}
        	ans=(ans*b)%p;
        	hashmap.Insert(ans,i);
    	}
    	
		ll t=Pow(b,m,p);
		ans=1;
    	for (int i=1;i<=m;i++){//gaint step
        	ans=(ans*t)%p;
        	if(hashmap.Find(ans)!=-1){
        		int t=i*m-hashmap.Find(ans)+cnt;
            	return(t%p+p)%p;
         	}
    	}
    	return -1;
	
}
int main(){
	while (scanf("%lld%lld%lld",&b,&p,&n)==3){  // b^x=n(mod p) 
		
		if (!b && !p && !n)
            break;
    	n%=p;
		ll tmp=bsgs(p,b,n);
		if(tmp==-1)printf("No Solution\n");
		else printf("%lld\n",tmp);

	}
}

---