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

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

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

定义:若aaa,ppp互质,则最小的正整数xxx使得ax1(modp)a^x \equiv 1 \pmod pax≡1(modp)即为aaa模ppp的阶。
性质xϕ(p)x | \phi(p)x∣ϕ(p)

原根

定义:若aaa模ppp的阶等于ϕ(p)\phi(p)ϕ(p),则称aaa为模ppp的一个原根。
当一个数aaa为模ppp的原根时,那么ai%pa^i \% pai%p的结果两两都不相同,且有2ap1,1ip1,iZ+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+
ppp有原根的充要条件为p=2,4,mn,2mnp = 2,4,m^n,2m^np=2,4,mn,2mn,其中mmm为奇质数,nnn是任意正整数(不会证)。
性质:当ppp有原根时,ppp的原根个数为ϕ(ϕ(p))\phi(\phi(p))ϕ(ϕ(p))。

指标

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

离散对数

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

BSGS算法

求解AxB(modP)A^x \equiv B \pmod PAx≡B(modP),其中AAA与PPP互质。
如果xxx有解,且xxx为所有解的最小正整数解,则0xϕ(P)0 \le x \le \phi(P)0≤x≤ϕ(P)。
根据欧拉定理,aϕ(p)1(modp)a^{ \phi(p) } \equiv 1 \pmod paϕ(p)≡1(modp),且a01(modp)a^0 \equiv 1 \pmod pa0≡1(modp),不难看出指数从000到ϕ(p)\phi(p)ϕ(p)为一个循环节(不一定最小哦(⊙o⊙))。
ϕ(p)\phi(p)ϕ(p)的最大值不超过p1p-1p−1,所以在000到p1p-1p−1一定可以找到一个可行的解。
我们设x=itjx = i * t - jx=i∗t−j,其中t=pt = \lceil \sqrt{p} \rceilt=⌈p​⌉(向下取整可能够不到ϕ(p)\phi(p)ϕ(p)),0jt10 \le j \le t - 10≤j≤t−1,0it0 \le i \le t0≤i≤t。
则方程变为AitjB(modP)A^{i * t - j} \equiv B \pmod PAi∗t−j≡B(modP)变为(At)iBAj(modP)(A^t)^i \equiv B * A^j \pmod P(At)i≡B∗Aj(modP),对于所有的0jt10 \le j \le t - 10≤j≤t−1,把BAj%pB*A^j \% pB∗Aj%p插入一个HASH表(可用map也可用邻接表)。
然后枚举iii的所有可能取值,计算出(At)i%p(A^t)^i \%p(At)i%p,在HASH表中查找是否存在对应的jjj,更新答案即可。时间复杂度为O(p)O(\sqrt p)O(p​)
还有特判就是如果A=0A=0A=0且B=0B=0B=0解为1,若A=0A=0A=0且BBB不等于0,无解。若BBB等于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算法

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

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++;
}

执行完代码后,原式变成DAxcntB(modP)D*A^{x-cnt} \equiv B \pmod {P}D∗Ax−cnt≡B(modP)我们继续设x=itj+cntx = i * t -j+cntx=i∗t−j+cnt,其中t=pt = \lceil \sqrt{p} \rceilt=⌈p​⌉,0jt10 \le j \le t - 10≤j≤t−1,0it0 \le i \le t0≤i≤t,用BSGS算法搞定。
因为i,j0i,j \ge 0i,j≥0,要求xcntx \ge cntx≥cnt,实际上存在x&lt;cntx \lt cntx<cnt的情况,直接套用会漏掉一些情况,对此我们要先进行一次O(log2P)O(log_2P)O(log2​P)的枚举(后面会讲到它的替代方法),从0到cnt1cnt - 1cnt−1枚举iii,直接验证Ai%P=BA^i \% P = BAi%P=B。
最后再说一些特判,如果在消去公因子的时候,出现B=DB=DB=D的情况,就已得出了答案,返回当时的cntcntcnt的值(这个操作与那个O(log2P)O(log_2P)O(log2​P)的枚举等价);如果B%PB \% PB%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;
}

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

解决第二个公式的方法(待填)

求解xAB(modP)x^A \equiv B \pmod PxA≡B(modP)

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

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

上一篇:费马小定理【证明】【复习】


下一篇:O(n)递推求乘法逆元