离散对数及其拓展 大步小步算法 BSGS

离散对数及其拓展

离散对数是在群ZpZ_{p}^{*}而言的,其中pp是素数。即在在群ZpZ_{p}^{*}内,aa是生成元,求关于xx的方程ax=ba^x=b的解,并将解记作x=logabx=log_{a}{b},离散对数指的就是这个logablog_{a}{b}.由于群ZpZ_{p}^{*}的阶是p1p-1,且是循环群,因为生成元的阶是p1p-1,因而模p1p-1相等的指数可以看做一样的数,x=logabx=log_{a}{b}的值不妨定为Zp1Z_{p-1}的元素,即模p1p-1的数,即0,1,2,3,,p20,1,2,3,\ldots,p-2.

以上说的是数学上对离散对数的定义,但是,搞算法竞赛人说的离散对数和这个有点差异,差异在于没有要求aa是生成元。

数学上的离散对数的定义,解一定存在且解模p1p-1意义下唯一,而我们平常说的离散对数由于没有要求aa是生成元,那么解可能不存在,解模p1p-1意义下不唯一,例如群Z7Z_7^*22不是生成元(因为23=12^3=1),于是以22为底的离散对数,模66意义下解不唯一,模33意义下解才唯一;而且当b=3,5,6b=3,5,6的时候没有解。

所以算法竞赛中说的离散对数一般是任意的aa,若有解求最小非负整数解,否则返回无解。

如何求离散对数的值?由于答案在Zp1Z_{p-1}内,最简单暴力的做法是枚举Zp1Z_{p-1}中的每个数,共p1p-1个数,复杂度O(p)O(p).

更高效的做法------Shank的大步小步算法

注意到ZpZ_{p}^{*}是个群,所以会有很多很好的性质,例如每个元素都有逆元。
我们选取一个数s,然后任意指数xx都可以表示成x=ks+r(0r<s)x=ks+r(0 \leq r < s).在群内,有------
ax=baks+r=b两边乘aksar=aksbar=(as)kb\begin{aligned} a^x &= b \\ a^{ks+r} &= b \quad\quad\quad\quad\quad \text{两边乘}a^{-ks}\\ a^r&=a^{-ks}b\\ a^r&= \left(a^{-s}\right)^{k}b \end{aligned}

每一个指数xx都可以和一个有序对(k,r)(k,r)建立一一对应关系,求xx就变成确定(k,r)(k,r)的问题了。注意上面的式子最后一行,左边ar(0r<k)a^r(0 \leq r < k)最多有ss个取值,右边aka^{-k}是固定的,而0kp2s0 \leq k \leq \frac{p-2}{s},因此右边最多有t=1+p2st=1+\lceil \frac{p-2}{s}\rceil个取值。

于是变成了右边tt个数在左边ss个数中有没有相等的元素.

这便是Shank的大步小步算法(Shank’s Baby-Step-Giant-Step
Algorithm)。左边的rr指数是每一加1,是小步;右边指数每次的变化是ss,是大步。

为了求最小的xx,具体地我们要从小到大枚举kk,计算出右边的值,查询左边ss个数中是否有相等的数。若有,符合条件最小的rrkk即构成最小的解xx;若无,则看下一个kk;若所有kk都看完了还没有,那么就说明无解。

为了加速查询,显然对左边的ss个数需要预先计算出来,并建立由值查询最小rr的索引"数组",如果直接用值做下标的话,空间需要O(p)O(p),但是总共才ss个值,所以应该用map建立"数组"或者hash(这个写麻烦一点)。

计算左边和右边所有取值的总的复杂度是O(s+t)O(s+t).插入及查询的时间复杂度使用map是O((s+t)lnr)O((s+t)\ln{r}),使用hash是O(s+t)O(s+t)。空间复杂度是索引数组大小O(s)O(s)。而t=1+p2st=1+\lceil \frac{p-2}{s}\rceil,取s=ps=\sqrt{p}附近的数,可以获得比较好的复杂度,如果使用map的话就是O(plnp)O(\sqrt{p}\ln{p}),如果是使用hash的话就是O(p)O(\sqrt{p}).
如果map会被卡就换成hash吧。

大步小步思路

略微思考便可以发现大步小步算法的精妙在于将原本需要的nn次计算的问题变成了左右各O(n)O(\sqrt{n})次的计算加O(n)O(\sqrt{n})次插入和查询问题。**这可以说是一种牺牲空间换时间的一种思路。**当我们去掉a,ba,b在这里代表的数含义的话,如果还可以有相同的变换,那么也可以如此降低复杂度。

另外,一点小优化,x=ks+rx=ks+r也可以换成x=ksr(k1,1rs)x=ks-r \quad(k \geq 1,1 \leq r \leq s)的形式,然后将式子变换成(as)k=bar(a^s)^k=ba^r.如此,就不用求逆元了。

会发现上面式子的变换只用到了群的性质,所以ZpZ_p^*换成ZnZ_n^*依旧可以如此做,只是阶不是由p1p-1变成了n1n-1,而是变成了ϕ(n)\phi(n),事实上p1p-1只是pp是质数时的ϕ(p)\phi(p)。当然,指数算到ϕ(n)1\phi(n)-1即可,当然,算到n2n-2也无妨,只是多算了一点而已,不影响结果。

当然,这就要求aa必须是ZnZ_n^*中的元素了,即(a,n)=1(a,n)=1,即互质。

拓展离散对数的思考与推导

离散对数的拓展,其实是想解axb(modn)a^x \equiv b \pmod{n}的方程中最小的自然数解x,其中aa是一般数,不保证与nn互质。因此必须思考aann不互质如何解的问题。

基本思路是,通过方程的等价变化,将暂时无法解决的问题化为可以解决的问题。

[如果只想看结论可以直接往下看下一小节的具体算法部分。]{.underline}

不妨考虑成qaxb(modn)qa^x \equiv b \pmod{n}形式的问题,首先将其视作关于axa^x的线性同余方程,线性同余方程有解等价于(q,n)b(q,n) \mid b.显然关于axa^x的线性同余方程有解是原方程有解的必要条件,此必要条件成立时,qaxb(modn)qa^x \equiv b \pmod{n}显然等价于qa(a,x)ax1b(a,n)(modn(a,n))\frac{qa}{(a,x)}a^{x-1} \equiv \frac{b}{(a,n)} \pmod{\frac{n}{(a,n)}}。新的方程和原方程形式上相同,但是模数变小了。如此,只要aa与模数不互质,就重复进行如此的变换,最终一定会终止于互质的情况(当然如果某一步等价变换的必要条件不满足则直接可以得出无解结论而提前终止了)。

假设最后终止于qcaxcbc(modnc)q_ca^{x-c} \equiv b_c \pmod{n_c},aann已经互质,则大步小步算法解之即可。

当指数xx在整个整数上域取值,这一系列的变换显然等价。

拓展离散对数的具体算法

解方程axb(modn)a^x \equiv b \pmod{n}的最小自然数解。

  1. 方程等价变换

    1. 将原方程通过等价变换,始终保持着qianibi(modni)q_{i}a^{n-i} \equiv b_i \pmod{n_i}的形式。

    2. q0=1,b0=b,n0=nq_0=1,b_0=b,n_0=n

    3. 迭代过程,qi+1=qia(a,ni),bi+1=bi(a,ni),ni+1=ni(a,ni)q_{i+1}=q_i\frac{a}{(a,n_i)},b_{i+1}=\frac{b_i}{(a,n_i)},n_{i+1}=\frac{n_i}{(a,n_i)},迭代的条件是bi+1=bi(a,ni)b_{i+1}=\frac{b_i}{(a,n_i)}是整数。

    4. 迭代终止条件:

      1. 不满足迭代条件返回无解

      2. (a,ni)=1(a,n_i)=1,终止迭代,假设终止时的ii取值是cc.

  2. 对迭代终止时的方程qcancbc(modnc)q_ca^{n-c} \equiv b_c \pmod{n_c}解关于anca^{n-c}的线性同余方程。

  3. 若线性同余方程无解,返回无解。

  4. 否则,假设解得ancB(modN)a^{n-c} \equiv B \pmod{N}

    1. (B,N)=1(B,N)=1,则a,Ba,B都在群ZNZ_{N}^{*}中,大步小步算法解axacB(modN)a^x \equiv a^cB \pmod{N}的最小自然数解x,终止。

    2. 否则,返回无解。

当然,最后先解线性同余方程,再进一步求x不是必要的。为了减少代码量,也可以直接再迭代终止形式的方程的基础上将aca^c调到右边,然后大步小步思想解决。事实上,大多数时候为了程序的简洁性都是如此做的。

Code

struct mod_sys{
	/*other mod_sys code*/
	// 这些方法和成员必须添加上
	set_mod,to_std and mod are needed 

	// 群Z_{p}^{*}下的离散对数log_a{b} 预设p=mod是素数
	// 使用大步小步算法
	// a^x = b (mod p)
	// 不要求a是生成元,但a mod p != 0,b mod p != 0
	// 由于时间复杂度 p^0.5*ln(p) 空间复杂度p^0.5
	// 故p^0.5肯定在int范围内,而且应该会更小一些
	// 预设p^2不爆ll ,否则乘法需要换成quick_mlt版本
	// 使用unordered_map作为查询数组
	// 返回是否有解,有解的情况下,x储存最小非负整数解
	bool log_p(ll a, ll b, ll &x) {
		a = to_std(a); b = to_std(b);
		unordered_map<ll,ll>val_to_r;
		ll s = (ll)ceil(sqrt((long long)mod));
		// x = ks-r r in [1,s] k in [1,(mod-2)/s+1]
		// (a^s)^k=b*a^r
		ll ar = 1,bar;
		for (ll r = 1; r <= s; ++r) {
			ar = (ar*a)%mod;
			bar = (b*ar)%mod;
			val_to_r[bar] = r; // 相同的k,查询应该返回最大的r,正序枚举r
		}
		// 循环结束,ar就是a^s
		ll &as = ar; 
		ll ask = 1; // (a^s)^k
		int t = (mod-2)/s+1;
		for (int k = 1; k <= t; ++k) {
			ask = (ask*as)%mod;
			auto it = val_to_r.find(ask);
			if (it != val_to_r.end()) {
				x = k*s-it->second;
				return true;
			}
		}
		return false;
	}

	// n=mod>=1,无其它特殊限制
	// (a,n)=1,(b,n)=1
	// a^x \equiv b (%n)
	// 返回是否有解
	// x存储最小自然数解
	bool log_n_easy(ll a,ll b, ll &x) {
		// x=ks-r in [0,phi_mod)
		// s \in [1,s], k in [1,(phi_mod-1)/s+1]
		ll phi_mod = phi(mod);
		a = to_std(a); b = to_std(b);
		unordered_map<ll,ll>val_to_r;
		ll s = (ll)ceil(sqrt((long long)phi_mod));
		ll ar = 1,bar; // a^r, b*a^r
		for (ll r = 1; r <= s; ++r) {
			ar = (ar*a)%mod;
			bar = (b*ar)%mod;
			val_to_r[bar] = r; // 相同的k,查询应该返回最大的r,正序枚举r
		}
		// 循环结束,ar就是a^s
		ll &as = ar; 
		ll ask = 1; // (a^s)^k
		// ask = bar
		int t = (phi_mod-1)/s+1;
		for (int k = 1; k <= t; ++k) {
			ask = (ask*as)%mod;
			auto it = val_to_r.find(ask);
			if (it != val_to_r.end()) {
				x = k*s-it->second;
				return true;
			}
		}
		return false;
	}

	// a^x = b (%mod) a,b,mod>=1没有特殊地限制
	// 返回是否有解,x存储最小自然数解
	bool log_n(ll a,ll b, ll &x) {
		ll q = 1,d,c=0;
		ll &n = mod;
		a = to_std(a); b = to_std(b);
		while (true) {
			if (q == b) return c;
			d = __gcd(a,n);
			if (1 == d) break;
			if (b%d) return false;
			b /= d; n /= d;
			q = a/d*(q%n)%n;
			a %= n;
			++c;
		}
		// qa^{x-c} = b (%n) <===>
		// 先解线性同余方程的步骤去掉,之后给log_n_easy添加一个参数q即可
		ll B,N;
		if (!linear_congruence_equation(q,b,n,B,N))
			return false;
		if (__gcd(B,N) != 1) return false;
		mod = N;
		bool t = log_n_easy(a,B,x);
		x += c;
		return t;
	}
};
posted @ 2019-09-11 21:03  我云知世就是力量  阅读(742)  评论(0)    收藏  举报