1. 引入
- Baby Step Giant Step算法(简称BSGS),用于求解形如 a x ≡ b ( m o d p ) a^x\equiv b(mod\ p) ax≡b(mod p)( a , b , p ∈ N a,b,p\in \mathbb{N} a,b,p∈N)的同余方程,即著名的离散对数问题。
- 本文分为 ( a , p ) = 1 (a,p)=1 (a,p)=1 和 ( a , p ) ≠ 1 (a,p)\neq 1 (a,p)̸=1 两种情况讨论。
2. 方程 a x ≡ b ( m o d p ) a^x\equiv b(mod\ p) ax≡b(mod p) 的解性
- 因为若 a x ≡ a x + n ( m o d p ) a^{x}\equiv a^{x+n}(mod\ p) ax≡ax+n(mod p),则 a x + i ≡ a x + n + i a^{x+i}\equiv a^{x+n+i} ax+i≡ax+n+i。
- 由抽屉原理以及同余的性质得, a x a^x ax 对 p p p 的模具有周期性,最大周期不超过 p p p。
- 由这个周期性可得,方程有解,等价于在 [ 0 , p ) [0,p) [0,p) 中有解。
- 由此我们只需要试图找到 [ 0 , p ) [0,p) [0,p) 中的最小自然数解,就可以得到方程是否有解,以及根据指数模的周期性得到方程的通解。
- 下面我们只讨论求解方程的最小自然数解。
3. 求解 a x ≡ b ( m o d p ) a^x\equiv b(mod\ p) ax≡b(mod p) ( a , p ) = 1 (a,p)=1 (a,p)=1
- Baby Step Giant Step它的名字就告诉我们,这个算法要将问题规模由大化小。
- 我们令 t = ⌈ p ⌉ t=\lceil \sqrt{p}\rceil t=⌈p⌉,那么我们所有 [ 0 , p ) [0,p) [0,p) 的自然数都可以包含在集合 { x ∣ x = i t − j , i ∈ [ 1 , t ] , j ∈ ( 0 , t ] } \{x|x=it-j,\ i\in [1,t],\ j \in (0,t]\} {x∣x=it−j, i∈[1,t], j∈(0,t]} 中。
- 我们只需要验证这个集合中是否存在解,即可验证原方程是否存在解。
- 我们考虑将方程表示为
a i t − j ≡ b ( m o d p ) a^{it-j}\equiv b(mod\ p) ait−j≡b(mod p)
- 根据 ( a , p ) = 1 (a,p)=1 (a,p)=1,我们有
a i t ≡ b a j ( m o d p ) a^{it}\equiv ba^j(mod\ p) ait≡baj(mod p)
- 我们发现等号两边的数我们都可以在 O ( p ) O(\sqrt p) O(p) 的时间内枚举出来。
- 我们只需要 O ( p ) O(\sqrt p) O(p) 枚举出所有 b a j m o d p ( j ∈ ( 0 , t ] ) ba^j\ mod\ p(\ j \in (0,t]) baj mod p( j∈(0,t]),然后把它们插入到哈希表当中,哈希表的对应位置储存这个值对应的最大的 j j j(因为要求最小自然数解,所以在 i i i 相同的时候要使 j j j 最大)。
- 然后 O ( p ) O(\sqrt p) O(p) 从小到大枚举我们需要的 a i t m o d p ( i ∈ [ 1 , t ] ) a^{it}\ mod\ p(\ i\in [1,t]) ait mod p( i∈[1,t]),并在哈希表中查询是否有相等的值,若存在则取最小的 i i i 及对应最大的 j j j,并将 i t − j it-j it−j 作为最小自然数解;否则若不存在则说明无自然数解。
inline int solve_BSGS(const int &a, const int &b, const int &p)
{
int t = ceil(sqrt(p));
std::map<int, int> hash;
//map实现hash表
hash.clear();
int tmp = b;
for (int i = 1; i <= t; ++i)
{
tmp = 1LL * tmp * a % p;
hash[tmp] = i;
}
//插入b*a^j
int pw = qpow(a, t, p);
tmp = pw;
for (int i = 1; i <= t; ++i)
{
if (hash.find(tmp) != hash.end())
return i * t - hash[tmp];
tmp = 1LL * tmp * pw % p;
}
//查询a^(it)
return -1; //返回无解
}
4. 求解 a x ≡ b ( m o d p ) a^x\equiv b(mod\ p) ax≡b(mod p) ( a , p ) ≠ 1 (a,p)\neq 1 (a,p)̸=1
- 对于这个方程,我们不能像上面那样求解的原因就是 a a a 在模 p p p 意义下不存在逆元,不能将 a i t − j a^{it-j} ait−j 表示为 a i t × a − j a^{it}\times a^{-j} ait×a−j。
- 那么我们从不定方程的角度分析这个同余方程。
- 这个方程等价于
a x + p y = b a^x+py=b ax+py=b
- 令 a 1 = ( a , p ) a_1=(a,p) a1=(a,p),令原方程化为
a x − 1 a a 1 + p a 1 y = b a 1 a^{x-1}\frac{a}{a_1}+\frac{p}{a_1}y=\frac{b}{a_1} ax−1a1a+a1py=a1b
- 若此时 ( a , p a 1 ) ≠ 1 (a,\frac{p}{a_1})\neq 1 (a,a1p)̸=1,那么我们接着化成
a x − 2 a 2 a 1 a 2 + p a 1 a 2 y = b a 1 a 2 a^{x-2}\frac{a^2}{a_1a_2}+\frac{p}{a_1a_2}y=\frac{b}{a_1a_2} ax−2a1a2a2+a1a2py=a1a2b
- 同理不断进行这样的操作,最后我们达到 ( a , p a 1 a 2 … a n ) = 1 (a,\frac{p}{a_1a_2\dots a_n})=1 (a,a1a2…anp)=1的目标,并将方程化为
a x − n a n a 1 a 2 … a n + p a 1 a 2 … a n y = b a 1 a 2 … a n a^{x-n}\frac{a^n}{a_1a_2\dots a_n}+\frac{p}{a_1a_2\dots a_n}y=\frac{b}{a_1a_2\dots a_n} ax−na1a2…anan+a1a2…anpy=a1a2…anb
- 然后记 a ′ = a a 1 a 2 … a n a'=\frac{a}{a_1a_2\dots a_n} a′=a1a2…ana, p ′ = p a 1 a 2 … a n p'=\frac{p}{a_1a_2\dots a_n} p′=a1a2…anp, b ′ = b a 1 a 2 … a n b'=\frac{b}{a_1a_2\dots a_n} b′=a1a2…anb,那么原不定方程可以化为同余方程
a x − n a ′ ≡ b ′ ( m o d p ′ ) a^{x-n}a'\equiv b'(mod\ p') ax−na′≡b′(mod p′)
- 显然 ( a ′ , p ′ ) = 1 (a',p')=1 (a′,p′)=1,因此我们可以写成
a x − n ≡ b ′ ( a ′ ) − 1 ( m o d p ′ ) a^{x-n}\equiv b'(a')^{-1}(mod\ p') ax−n≡b′(a′)−1(mod p′)
-
然后就可以用互质的方法解决了。
-
可以发现,每次在 p p p 中除去一个最大公约数,每次都会有至少同一个质因子的次数减少 1 1 1,那么在
int
范围内, n n n 最多只会取到 30 30 30。 -
以上两个情况的时间复杂度均为 O ( p log 2 p ) O(\sqrt{p}\log_2{\sqrt p}) O(plog2p),因为用map实现哈希,可以做到更优秀。
-
有一些值得注意的地方:
- 在我们令 b b b 除以一个最大公约数 d d d 时,若 d ∤ b d\nmid b d∤b,结合求解二元不定方程的知识,我们判定方程无自然数解。
- 我们在不互质情况的化简操作中,已经假定了 x ≥ n x\ge n x≥n,所以方程才能写成那样的互质形式。对于 x ≤ n x\le n x≤n 的情况,我们应当提前枚举判断。
- 注意在求解方程之前将 a , b a,b a,b 对 p p p 取模,注意 a a a 取模后为 0 0 0 的情况。
- 对于对时间限制要求较为紧的题目,应当使用更为优秀的哈希表实现方式。
-
模板题:SP3105 MOD - Power Modulo Inverted
-
代码:
//map水过
#include <map>
#include <cmath>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
int a, p, b;
inline int qpow(int b, int p, const int &mod)
{
int res = 1;
for (; p; p >>= 1, b = 1LL * b * b % mod)
if (p & 1)
res = 1LL * res * b % mod;
return res;
}
inline int ex_gcd(const int &a, const int &b, int &x, int &y)
{
if (!b)
return x = 1, y = 0, a;
int res = ex_gcd(b, a % b, y, x);
return y -= a / b * x, res;
}
inline int solve_equ(const int &a, const int &b, const int &c)
{
int x, y;
int d = ex_gcd(a, b, x, y);
if (c % d != 0) return -1;
int mod = b / d;
return (1LL * c / d * x % mod + mod) % mod;
}
inline int solve_BSGS(const int &a, const int &b, const int &p)
{
int t = ceil(sqrt(p));
std::map<int, int> hash;
hash.clear();
int tmp = b;
for (int i = 1; i <= t; ++i)
{
tmp = 1LL * tmp * a % p;
hash[tmp] = i;
}
int pw = qpow(a, t, p);
tmp = pw;
for (int i = 1; i <= t; ++i)
{
if (hash.find(tmp) != hash.end())
return i * t - hash[tmp];
tmp = 1LL * tmp * pw % p;
}
return -1;
}
inline bool check()
{
int k = 1 % p;
for (int i = 0; i <= 40; ++i)
{
if (k == b)
return printf("%d\n", i), true;
k = 1LL * k * a % p;
}
if (!a)
return puts("No Solution"), true;
return false;
}
int main()
{
while (scanf("%d%d%d", &a, &p, &b), a || p || b)
{
a %= p, b %= p;
if (check())
continue;
int d;
int ap = 1, n = 0;
bool flg = false;
while ((d = std::__gcd(a, p)) != 1)
{
++n;
ap = 1LL * ap * (a / d) % p;
p /= d;
if (b % d)
{
flg = true;
break;
}
b /= d;
}
if (flg)
puts("No Solution");
else
{
int res = solve_BSGS(a, 1LL * b * solve_equ(ap, p, 1) % p, p);
if (res == -1)
puts("No Solution");
else
printf("%d\n", res + n);
}
}
return 0;
}