同余3：解高次同余方程的BSGS算法及其扩展学习笔记

• 阶
• 原根
• 指标
• 离散对数
• BSGS算法
• 扩展BSGS算法
• 解决第二个公式的方法（待填）

前言：在前两篇博客，我提到了解决单个线性同余方程的方法，以及解决线性同余方程组的方法，但是，当单个同余方程变成$A^x \equiv B \pmod P$Ax≡B(modP)或$x^A \equiv B \pmod P$xA≡B(modP)时，蒟蒻的我又得自闭了，为了不再自闭，就不得不提到BSGS算法 。BSGS算法，原名Baby Step,Giant Step算法，不过，有人称之为北上广深算法，或者拔山盖世算法。经过本蒟蒻的鉴定，这个算法需要用到一些欧拉定理的知识，所以，我们不妨先跑一下题#手动滑稽# 。
本文涉及到的欧拉定理：当$a$a与$p$p互质时，有$a^{ \phi(p) } \equiv 1 \pmod p$aϕ(p)≡1(modp)

阶

定义：若$a$a，$p$p互质，则最小的正整数$x$x使得$a^x \equiv 1 \pmod p$ax≡1(modp)即为$a$a模$p$p的阶。
性质：$x | \phi(p)$x∣ϕ(p)

原根

定义：若$a$a模$p$p的阶等于$\phi(p)$ϕ(p)，则称$a$a为模$p$p的一个原根。
当一个数$a$a为模$p$p的原根时，那么$a^i \% p$ai%p的结果两两都不相同，且有$2 \le a \le p - 1,1 \le i \le p - 1,i \in Z^+$2≤a≤p−1,1≤i≤p−1,i∈Z+
模$p$p有原根的充要条件为$p = 2,4,m^n,2m^n$p=2,4,mn,2mn，其中$m$m为奇质数，$n$n是任意正整数（不会证）。
性质：当$p$p有原根时，$p$p的原根个数为$\phi(\phi(p))$ϕ(ϕ(p))。

指标

如果$a$a是$p$p的原根，对于任意一个整数$x$x小于$p$p，都可以找到一个$l(x)$l(x)小于等于$\phi(m)$ϕ(m)，使得$a^{l(x)} \equiv x \pmod p$al(x)≡x(modp)，则$l(x)$l(x)称为$x$x的指标。
法则1：$l(ab) \equiv l(a) + l(b) \pmod {\phi(p)}$l(ab)≡l(a)+l(b)(modϕ(p))
法则2：$l(a^k) \equiv k*l(a) \pmod {\phi(p)}$l(ak)≡k∗l(a)(modϕ(p))

离散对数

一种在整数中基于同余运算和原根的对数运算。当模$p$p有原根时，设$a$a为模$p$p的一个原根，则当$x \equiv a^p \pmod m$x≡ap(modm)时，$log_ax \equiv p \pmod {\phi(m)}$loga​x≡p(modϕ(m))，此处的$log_ax$loga​x为以$a$a为底，模$\phi(m)$ϕ(m)的离散对数值。

BSGS算法

求解$A^x \equiv B \pmod P$Ax≡B(modP)，其中$A$A与$P$P互质。
如果$x$x有解，且$x$x为所有解的最小正整数解，则$0 \le x \le \phi(P)$0≤x≤ϕ(P)。
根据欧拉定理，$a^{ \phi(p) } \equiv 1 \pmod p$aϕ(p)≡1(modp)，且$a^0 \equiv 1 \pmod p$a0≡1(modp)，不难看出指数从$0$0到$\phi(p)$ϕ(p)为一个循环节（不一定最小哦(⊙o⊙)）。
而$\phi(p)$ϕ(p)的最大值不超过$p-1$p−1，所以在$0$0到$p-1$p−1一定可以找到一个可行的解。
我们设$x = i * t - j$x=i∗t−j，其中$t = \lceil \sqrt{p} \rceil$t=⌈p​⌉（向下取整可能够不到$\phi(p)$ϕ(p)），$0 \le j \le t - 1$0≤j≤t−1，$0 \le i \le t$0≤i≤t。
则方程变为$A^{i * t - j} \equiv B \pmod P$Ai∗t−j≡B(modP)变为$(A^t)^i \equiv B * A^j \pmod P$(At)i≡B∗Aj(modP)，对于所有的$0 \le j \le t - 1$0≤j≤t−1，把$B*A^j \% p$B∗Aj%p插入一个HASH表（可用map也可用邻接表）。
然后枚举$i$i的所有可能取值，计算出$(A^t)^i \%p$(At)i%p，在HASH表中查找是否存在对应的$j$j，更新答案即可。时间复杂度为$O(\sqrt p)$O(p​)
还有特判就是如果$A=0$A=0且$B=0$B=0解为1，若$A=0$A=0且$B$B不等于0，无解。若$B$B等于1，直接返回零（不知为何少了这个特判会错）。
板子http://poj.org/problem?id=2417

#include <iostream>
#include <cstdio>
#include <cstring>
#include <map>
#include <cmath>
#define ll long long
using namespace std;
inline int power( int a , int b , int p ) {
	int res = 1;
	a %= p; 
	while (b) {
		if ( b & 1 ) {
			res = ( (ll)res * (ll)a ) % p;
		}
		a = ( (ll)a * (ll)a ) % p;
		b >>= 1;
	}
	return res;
}
inline int bsgs( int a , int b , int p ) {//若无解则返回-1 
	if ( b == 1 ) {
		return 0;//特判 
	}
	if ( !a ) {
		return b == 0 ? 1 : -1;//特判
	}
	map< int , int > hash;
	hash.clear();
	b %= p;
	int t = (int)(sqrt(p * 1.0)) + 1;
	for ( int j = 0 ; j <= t - 1 ; ++j ) {
		int val = ( (ll)b * (ll)power( a , j , p ) ) % p;
		hash[val] = j;
	}
	a = power( a , t , p );
	for ( int i = 0 ; i <= t ; ++i ) {
		int val = power( a , i , p );
		if ( hash.find(val) == hash.end() ) {
			continue;
		}
		int temp = hash[val];
		if ( i * t - temp >= 0 ) {
			return i * t - temp;
		}
	}
	return -1;
}
int main () {
	int a , b , p;
	while ( scanf("%d%d%d",&p,&a,&b) != EOF ) {
		int tem = bsgs( a , b , p );
		if ( tem == -1 ) {
			printf("no solution\n");
		} else {
			printf("%d\n",tem);
		}
	}
	return 0;
}

map在poj上会T飞，考虑邻接表+真正的hash。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#define ll long long
#define mod 100007
using namespace std;
struct hashtable{
	int bcnt;
	int head[mod];
	struct line{
		int next;
		int to;
		int val;
	}str[10 * mod];
	inline void insert( int from , int to , int val ) {
		str[++bcnt].next = head[from];
		head[from] = bcnt;
		str[bcnt].to = to;
		str[bcnt].val = val;
		return;
	}
	inline void clear() {
		bcnt = 0;
		memset( head , 0 , sizeof(head) );
	}
	inline void doit( int x , int y ) {
		int te = x % mod;
		insert( te , y , x );
	}
	inline int query( int x ) {
		int te = x % mod;
		for ( int i = head[te] ; i ; i = str[i].next ) {
			int sn = str[i].to;
			int va = str[i].val;
			if ( va == x ) {
				return sn;
			}
		}
		return -1;
	}
}hash;
inline int power( int a , int b , int p ) {
	int res = 1;
	a %= p; 
	while (b) {
		if ( b & 1 ) {
			res = ( (ll)res * (ll)a ) % p;
		}
		a = ( (ll)a * (ll)a ) % p;
		b >>= 1;
	}
	return res;
}
inline int bsgs( int a , int b , int p ) {//若无解则返回-1 
	if ( b == 1 ) {
		return 0;//特判 
	}
	if ( !a ) {
		return b == 0 ? 1 : -1;//特判
	}
	hash.clear();
	b %= p;
	int t = (int)(sqrt(p * 1.0)) + 1;
	for ( int j = 0 ; j <= t - 1 ; ++j ) {
		int val = ( (ll)b * (ll)power( a , j , p ) ) % p;
		hash.doit( val , j );
	}
	a = power( a , t , p );
	for ( int i = 0 ; i <= t ; ++i ) {
		int val = power( a , i , p );
		int temp = hash.query(val);
		if ( temp == -1 ) {
			continue;
		}
		if ( i * t - temp >= 0 ) {
			return i * t - temp;
		}
	}
	return -1;
}
int main () {
	int a , b , p;
	while ( scanf("%d%d%d",&p,&a,&b) != EOF ) {
		int tem = bsgs( a , b , p );
		if ( tem == -1 ) {
			printf("no solution\n");
		} else {
			printf("%d\n",tem);
		}
	}
	return 0;
}

扩展BSGS算法

还是求解$A^x \equiv B \pmod P$Ax≡B(modP)，其中$A$A与$P$P不互质。
大致做法跟BSGS算法差不多，只是多了一个消除因子的操作，不难想到：
设$A=a*d$A=a∗d，$B=b*d$B=b∗d，$P=p*d$P=p∗d。
对于原式$(a*d)^x \equiv b*d \pmod {p*d}$(a∗d)x≡b∗d(modp∗d)根据同余的一条性质，原式等价于$a*(a*d)^{x-1} \equiv b \pmod p$a∗(a∗d)x−1≡b(modp)因而，开始时我们不断消除$A$A，$P$P的公因子，为了方便，采取将$P$P，$B$B不断除以公因子，对于$A^x$Ax记录被消去公因子的$A$A的个数，如果$B$B不能整除这些公因子(其实还有一些特判)，那么问题无解。
其实上述操作等价于消除$A^x$Ax，$P$P的公因子，如果觉得我的表述晦涩难懂，可以参考以下代码。
我以后的表述会用上以下代码中的变量。

ll D = 1 , cnt = 0;
while ( gcd( A , P ) != 1 ) {
	if ( B % gcd( A , P ) != 0 ) {
		printf("No Solution");
		return 0;
	}
	B /= gcd( A , C );
	P /= gcd( A , C );
	D *= A / gcd( A , C );
	cnt++;
}

执行完代码后，原式变成$D*A^{x-cnt} \equiv B \pmod {P}$D∗Ax−cnt≡B(modP)我们继续设$x = i * t -j+cnt$x=i∗t−j+cnt，其中$t = \lceil \sqrt{p} \rceil$t=⌈p​⌉，$0 \le j \le t - 1$0≤j≤t−1，$0 \le i \le t$0≤i≤t，用BSGS算法搞定。
因为$i,j \ge 0$i,j≥0，要求$x \ge cnt$x≥cnt，实际上存在$x \lt cnt$x<cnt的情况，直接套用会漏掉一些情况，对此我们要先进行一次$O(log_2P)$O(log2​P)的枚举（后面会讲到它的替代方法），从0到$cnt - 1$cnt−1枚举$i$i，直接验证$A^i \% P = B$Ai%P=B。
最后再说一些特判，如果在消去公因子的时候，出现$B=D$B=D的情况，就已得出了答案，返回当时的$cnt$cnt的值(这个操作与那个$O(log_2P)$O(log2​P)的枚举等价)；如果$B \% P$B%P等于1，直接返回0，其余特判跟BSGS算法相同。
这里我就不用map了，直接上邻接表。
板子http://poj.org/problem?id=3243

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#define ll long long
#define mod 100007
using namespace std;
struct hashtable{
	int bcnt;
	int head[mod];
	struct line{
		int next;
		int to;
		int val;
	}str[10 * mod];
	inline void insert( int from , int to , int val ) {
		str[++bcnt].next = head[from];
		head[from] = bcnt;
		str[bcnt].to = to;
		str[bcnt].val = val;
		return;
	}
	inline void clear() {
		bcnt = 0;
		memset( head , 0 , sizeof(head) );
	}
	inline void doit( int x , int y ) {
		int te = x % mod;
		insert( te , y , x );
	}
	inline int query( int x ) {
		int te = x % mod;
		for ( int i = head[te] ; i ; i = str[i].next ) {
			int sn = str[i].to;
			int va = str[i].val;
			if ( va == x ) {
				return sn;
			}
		}
		return -1;
	}
}hash;
inline int gcd( int a , int b ) {
	return b == 0 ? a : gcd( b , a % b );
}
inline int power( int a , int b , int p ) {
	int res = 1;
	a %= p; 
	while (b) {
		if ( b & 1 ) {
			res = ( (ll)res * (ll)a ) % p;
		}
		a = ( (ll)a * (ll)a ) % p;
		b >>= 1;
	}
	return res;
}
void c() {
	printf("yes\n");
}
inline int exbsgs( int a , int b , int p ) {//若无解则返回-1 
	a %= p;
	b %= p;
	if ( b == 1 ) {
		return 0;//特判 
	}
	if ( !a ) {
		return b == 0 ? 1 : -1;//特判
	}
	int g , cnt = 0;
	ll d = 1;
	while ( gcd( a , p ) > 1 ) {
		g = gcd( a , p );
		if ( b % g ) {
			return -1;
		}
		b /= g;
		p /= g;
		++cnt;
		d = d * ( a / g ) % p;
		if ( b == d ) {
			return cnt;
		}
	}
	hash.clear();
	int t = (int)(sqrt(p * 1.0)) + 1;
	for ( int j = 0 ; j <= t - 1 ; ++j ) {
		int val = ( (ll)b * (ll)power( a , j , p ) ) % p;
		hash.doit( val , j );
	}
	a = power( a , t , p );
	for ( int i = 0 ; i <= t ; ++i ) {
		int val = power( a , i , p );
		val = ((ll)val * (ll)d) % p;//不要落了d
		int temp = hash.query(val);
		if ( temp == -1 ) {
			continue;
		}
		if ( i * t - temp + cnt >= 0 ) {
			return i * t - temp + cnt;
		}
	}
	return -1;
}
int main () {
	int a , b , p;
	while ( scanf("%d%d%d",&a,&p,&b) != EOF && a && b && p ) {
		int tem = exbsgs( a , b , p );
		if ( tem == -1 ) {
			printf("No Solution\n");
		} else {
			printf("%d\n",tem);
		}
	}
	return 0;
}

讲了这么多，我们终于得到了$A^x \equiv B \pmod P$Ax≡B(modP)的解法，那么第二个公式怎么办(⊙﹏⊙)？难道又要自闭吗 ？

解决第二个公式的方法（待填）

求解$x^A \equiv B \pmod P$xA≡B(modP)

这需要用到一些数论知识，等数论书到了再填坑吧。

本文参考资料：《算法竞赛进阶指南》
《信息学竞赛之数学一本通》
百度百科