【学习笔记】exBSGS

众所周知,朴素的 BSGS 并不可以解决 $(a,p) >1$ 的情况。然后 exBSGS 大概就是魔改一下 BSGS。

Preface

严格来讲,这东西不应该被算到一个模板里面。因为在我看来模板是人构造出来的,但是这个算法应该是一个解决问题的process…更像是在解一道数学题,如果 BSGS 是定理的话,exBSGS 更像是一个不断转化的过程233(手动@lxa并且溜

Algorithm~Process

今天才发现原来$\rm{BSGS}$有两种写法…(PS on 2020.4.16 : 明明有好多好多种写法)

其实本质上,当 $p$ 不为素数时,我们无法进行朴素BSGS 的原因是我们的欧拉定理 $a^{\varphi(p)} \equiv b(\bmod p)$ 只能处理 $(a,p)=1$ 的情况。那么朴素的 $\rm{BSGS}$ 关键在于,可以保证最小解是有界的—— $x$ 一定在 $[1,\varphi(p)]$ 中。所以最后 $\rm BSGS$ 的复杂度才会是 $\Theta(\sqrt{\varphi(p)})$ 的——比如说比较常见的 $p$ 是素数的情况下,时间复杂度为 $\Theta(\sqrt p)$ 。

那么也就是说,我们只需要进行一些操作,保证$(a,p)=1 $即可$^{[1]}$。

对于同余式 $a^x\equiv b \pmod p$ 而言,先假定 $(a,p)>1 $。而此时如果有 $((a,p), b)=1$,那么说明此式只有可能在 $x=0,b=1$ 的时候有解——这个结论是平凡的。因为假设我们把它展开成 $a\cdot a^{x-1} +kp=b $ 的形式,必须要有$(a,p) | b$ 的情况下,才能保证 $a^{x-1}$ 和 $k$ 都是整数。

那么对于 $(a,p)>1$ 且 $(a,p) | b $,将原式变成

的样子,如果此时 $(a^{x-1},\frac{p}{(a,p)})=1$ 的话,就直接解

这个方程即可。否则我们继续分解直至 $(p’,a)=1$。

那么此时有个问题需要注意,就是如果在解这个方程时,出现了

的情况,那我们需要特判并 return -1 ,因为此时 $(a^{x-2},\frac{p}{(a,p)})$ 不存在整数解。

另一种情况,如果我们出现了

的情况,也需要特判并输出此 $k$(此时同余式左边是 $a^{x-k}$,因为 $a^{x-k}\equiv1 \pmod p$ 所以直接输出 $k$ ),不过也有可能不需要,完全看你写的 BSGS 能不能判断 $x=0$ 的情况…一般情况下不能。

此时由于 $\boldsymbol{p}$ 不再是素数,所以不能用费马小定理,需要我们用 exgcd 的方法求逆元,包括但不限于 $\frac{b}{(a,p)}$ 的逆元和 $a^{-im}$

以下是完整版代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <map>
#include <cmath>
#include <cstdio>
#include <iostream>
#include <unordered_map>

#define ll long long

using namespace std ;

unordered_map<ll, int> H ;
int N, M, P, ans ; // N ^x = M (mod P)

inline ll gcd(ll a, ll b){
if (!b) return a ;
return gcd(b, a % b) ;
}
inline ll expow(ll a, ll b, ll mod){
ll res = 1 ;
while (b) res = ((b & 1)?res * a % mod : res), a = a * a % mod, b >>= 1 ;
return res ;
}
inline ll exgcd(ll &x, ll &y, ll a, ll b){
if (!b){ x = 1, y = 0 ; return a ; }
ll t = exgcd(y, x, b, a % b) ; y -= x * (a / b) ; return t ;
}
inline ll BSGS(ll a, ll b, ll mod, ll qaq){
H.clear() ; ll Q, p = ceil(sqrt(mod)), x, y ;
exgcd(x, y, qaq, mod), b = (b * x % mod + mod) % mod,
Q = expow(a, p, mod), exgcd(x, y, Q, mod), Q = (x % mod + mod) % mod ;
for (ll i = 1, j = 0 ; j <= p ; ++ j, i = i * a % mod) if (!H.count(i)) H[i] = j ;
for (ll i = b, j = 0 ; j <= p ; ++ j, i = i * Q % mod) if (H[i]) return j * p + H[i] ; return -1 ;
}
inline ll exBSGS(){
ll qaq = 1 ;
ll k = 0, qwq = 1 ;
if (M == 1) return 0 ;
while ((qwq = gcd(N, P)) > 1){
if (M % qwq) return -1 ;
++ k, M /= qwq, P /= qwq, qaq = qaq * (N / qwq) % P ;
if (qaq == M) return k ;
}
return (qwq = BSGS(N, M, P, qaq)) == -1 ? -1 : qwq + k ;
}
int main(){
while(cin >> N){
scanf("%d%d", &P, &M);
if (!N && !M && !P) return 0 ;
N %= P, M %= P, ans = exBSGS() ;
if (ans < 0) puts("No Solution") ; else cout << ans << '\n' ;
}
}

然后 2020.4.16 我又重写了一遍,似乎是比原来快了…但是还是比 yyb 慢得多啊。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <unordered_map>

using namespace std ;

#define il inline

#define fuck_out puts("No Solution")

typedef long long ll ;

int T ;
int ans ;

namespace exBSGS{
unordered_map <int, int> f ;
il int expow(int a, int b, int m){
int ret = 1 ;
if (a > m) a -= m ;
while (b){
if (b & 1)
ret = 1ll * ret * a % m ;
a = 1ll * a * a % m ; b >>= 1 ;
}
return ret ;
}
il bool chk_ans(int &a, int b, int m){
a %= m ;
if (b >= m) return 1 ;
if (!a && b) return 1 ;
return 0 ;
}
int gcd(int x, int y){
if (!y) return x ;
return gcd(y, x % y);
}
il int phi(int x){
int y = sqrt(x), ret = x ;
for (int i = 2 ; i <= y ; ++ i){
if (x % i == 0)
ret = 1ll * ret * (i - 1) / i ;
while (x % i == 0) x /= i ;
}
if (x > 1)
ret = 1ll * ret * (x - 1) / x ;
return ret ;
}
int exgcd(int& x, int& y, int a, int b){
if (!b){
x = 1 ;
y = 0 ;
return a ;
}
int ret = exgcd(y, x, b, a % b) ;
y = y - (ll)x * (a / b) ; return ret ;
/*
exgcd(x, y, b, a % b) ;
int t = y ; y = x ;
x = x - (ll)y * a / b ;
*/
}
il int bsgs(int a, int b, int m){
int p = sqrt(m) + 1 ;
int bl ; int t = - p ;
while (t < 0) t += (m - 1) ;
bl = expow(a, t, m) ; f.clear() ;
if (chk_ans(a, b, m)) return -1 ;
for (int i = 0, j = 1 ; i < p ; ++ i){
if (!f.count(j)) f[j] = i ; j = (ll)j * a % m ;
}
for (int i = 0, j = b ; i <= p ; ++ i){
if (f.count(j)) return f[j] + i * p ; j = (ll)j * bl % m ;
}
return -1 ;
}
il int bsgs(int a, int b1, int m, int b2){
f.clear() ;
int x, y, p, bl ;
exgcd(x, y, b2, m) ;
b1 = ((ll)b1 * x % m + m) % m ;
if (chk_ans(a, b1, m)) return -1 ;
exgcd(x, y, a, m) ; (x %= m) += m ;
p = ceil(sqrt(phi(m))) ; bl = expow(x, p, m) ;
for (int i = 0, j = 1 ; i < p ; ++ i){
if (!f.count(j)) f[j] = i ; j = (ll)j * a % m ;
}
for (int i = 0, j = b1 ; i <= m / p ; ++ i){
if (f.count(j)) return f[j] + i * p ; j = (ll)j * bl % m ;
}
return -1 ;
}
il int exbsgs(int a, int b, int m){
if (m <= 1) return -1 ;
int g, ans = 0, res, inv = 1 ;
while ((g = gcd(a, m)) > 1){
if (b % g != 0) return -1 ;
b /= g ; m /= g ; ans += 1 ;
inv = 1ll * inv * (a / g) % m ;
if (inv == m) return ans ;
}
res = bsgs(a, b, m, inv) ;
if (res < 0) return -1 ; return res + ans ;
}
}
using namespace exBSGS ;

int main(){
int n, m, k ;
while (1){
scanf("%d%d%d", &n, &k, &m) ;
if (!(n + k + m)) return 0 ; ans = exbsgs(n, m, k) ;
if (ans < 0) fuck_out ; else printf("%d\n", ans) ;
}
return 0 ;
}

Afterword

今天才知道原来 BSGS 有两种写法qaq

$zyf2000$ 好像和我写的 BSGS 对“大步”和“小步”的定义不是很一样…但其实正着做或者反着做复杂度都是对的。

于是最后还是自己$\rm{yy}$的233

不过整体来看其实不是很难。因为本质上 exBSGS 就是 BSGS 和「同余」那一节里面的一个小定理结合了一下而已。

$\rm{Reference}$