BSGS算法学习笔记
从这里开始
离散对数和BSGS算法
设$x$是最小的非负整数使得$a^{x}\equiv b\ \ \ \pmod{m}$，则$x$是$b$以$a$为底的离散对数，记为$x = ind_{a}b$。
假如给定$a, b, m$，考虑如何求$x$，或者输出无解，先考虑$(a, m) = 1$的情况。
定理1（欧拉定理） 若$(a, m) = 1$，则$a^{\varphi(m)}\equiv 1 \pmod{m}$。
证明这里就不给出，因为在百度上随便搜一搜就能找到。
不过，这个定理告诉我们，在$(a, m) = 1$的情况下，若存在答案，则答案不会超过$\varphi(m) - 1$。
考虑$a^{x} \equiv b \pmod{m}$，通过一些操作可以得到：
$a^{x - k} \equiv a^{-k}b \pmod{m}$
因此可以选取正整数$c$，将$x$表示为$ic + j$的形式，然后有：
$a^{ic} \equiv a^{-j}b \pmod{m}$
考虑预处理$a^{-j}b$，以它的值为键，最小的$j$为值存入Hash表或者Map中。
这样有什么用呢？你可以快速枚举$a^{ic}$，然后你将这个值在Hash表中查一查对应的最小的$j$，如果查到就可以得到答案了。
Code
1 /** 2 * poj 3 * Problem#2417 4 * Accepted 5 * Time: 16ms 6 * Memory: 1372l 7 */ 8 #include <iostream> 9 #include <cstring> 10 #include <cstdio> 11 #include <cmath> 12 using namespace std; 13 typedef bool boolean; 14 15 int p, x, a; 16 17 typedef class HashMap { 18 private: 19 static const int M = 46666; 20 public: 21 int ce; 22 int h[M], key[M], val[M], next[M]; 23 24 HashMap():ce(-1) { } 25 26 void insert(int k, int v) { 27 int ha = k % M; 28 for (int i = h[ha]; ~i; i = next[i]) 29 if (key[i] == k) { 30 val[i] = v; 31 return; 32 } 33 ++ce, key[ce] = k, val[ce] = v, next[ce] = h[ha]; 34 h[ha] = ce; 35 } 36 37 int operator [] (int k) { 38 int ha = k % M; 39 for (int i = h[ha]; ~i; i = next[i]) 40 if (key[i] == k) 41 return val[i]; 42 return -1; 43 } 44 45 void clear() { 46 ce = -1; 47 memset(h, -1, sizeof(h)); 48 } 49 }HashMap; 50 51 int qpow(int a, int pos) { 52 int pa = a, rt = 1; 53 for (; pos; pos >>= 1, pa = pa * 1ll * pa % p) 54 if (pos & 1) 55 rt = rt * 1ll * pa % p; 56 return rt; 57 } 58 59 void exgcd(int a, int b, int& d, int &x, int &y) { 60 if (!b) 61 d = a, x = 1, y = 0; 62 else { 63 exgcd(b, a % b, d, y, x); 64 y -= (a / b) * x; 65 } 66 } 67 68 int inv(int a, int n) { 69 int d, x, y; 70 exgcd(a, n, d, x, y); 71 return (x < 0) ? (x + n) : (x); 72 } 73 74 inline boolean init() { 75 return ~scanf("%d%d%d", &p, &x, &a); 76 } 77 78 int cs; 79 HashMap mp; 80 inline int ind() { 81 mp.clear(); 82 cs = sqrt(p - 1 + 0.5); 83 if (cs == 0) cs++; 84 int ainv = inv(x, p), iap = a * 1ll * qpow(ainv, cs - 1) % p; 85 for (int i = cs - 1; ~i; i--, iap = iap * 1ll * x % p) 86 mp.insert(iap, i); 87 int cp = qpow(x, cs), pw = 1; 88 for (int i = 0; i < p; i += cs, pw = pw * 1ll * cp % p) 89 if (~mp[pw]) 90 return mp[pw] + i; 91 return -1; 92 } 93 94 inline void solve() { 95 int res = ind(); 96 if (res == -1) 97 puts("no solution"); 98 else 99 printf("%d\n", res); 100 } 101 102 int main() { 103 while (init()) 104 solve(); 105 return 0; 106 }
扩展BSGS算法
使刚刚的问题更一般，去掉$(a, m) = 1$的条件。
此时逆元不一定存在，所以不能用上述的BSGS算法来做。
考虑去掉它们的公约数$d$。
得到
$a^{x - 1}\cdot\frac{a}{d} \equiv \frac{b}{d} \pmod{\frac{m}{d}}$
在这步中，如果$b \nmid d$，那么显然无解。
否则我可以令$x' = x - 1,k' = \frac{a}{d}, b'=\frac{b}{d}, m' = \frac{m}{d}$进行换元得到：
$k'a^{x'} \equiv b' \pmod{m'}$
如果$(a, m') = 1$，那么直接BSGS解这个方程，然后带回去算原先的$x$，否则可以继续计算$a$和$m'$的最大公约数，继续除掉它，直到$(a, m) = 1$，然后BSGS解方程。
因为最大公约数不为1，每次至少除以2，1和任何数互质，因此总共除的次数不会超过$\log_{2}m$。
但是这么做存在一个问题，假如除的次数为$k$，那么它会忽略大于等于0小于$k$的解，因此，除的时候判一下即可。
Code
1 /** 2 * poj 3 * Problem#3243 4 * Accepted 5 * Time: 47ms 6 * Memory: 1248k 7 */ 8 #include <iostream> 9 #include <cstring> 10 #include <cstdio> 11 #include <cmath> 12 #ifndef WIN32 13 #define Auto "%lld" 14 #else 15 #define Auto "%I64d" 16 #endif 17 using namespace std; 18 typedef bool boolean; 19 20 int x, a, m; 21 22 typedef class HashMap { 23 private: 24 static const int M = 46666; 25 public: 26 int ce; 27 int h[M], key[M], val[M], next[M]; 28 29 HashMap():ce(-1) { } 30 31 void insert(int k, int v) { 32 int ha = k % M; 33 for (int i = h[ha]; ~i; i = next[i]) 34 if (key[i] == k) { 35 val[i] = v; 36 return; 37 } 38 ++ce, key[ce] = k, val[ce] = v, next[ce] = h[ha]; 39 h[ha] = ce; 40 } 41 42 int operator [] (int k) { 43 int ha = k % M; 44 for (int i = h[ha]; ~i; i = next[i]) 45 if (key[i] == k) 46 return val[i]; 47 return -1; 48 } 49 50 void clear() { 51 ce = -1; 52 memset(h, -1, sizeof(h)); 53 } 54 }HashMap; 55 56 int qpow(int a, int pos, int m) { 57 int pa = a, rt = 1; 58 for (; pos; pos >>= 1, pa = pa * 1ll * pa % m) 59 if (pos & 1) 60 rt = rt * 1ll * pa % m; 61 return rt; 62 } 63 64 int gcd (int a, int b) { 65 return (b) ? (gcd(b, a % b)) : (a); 66 } 67 68 void exgcd(int a, int b, int& d, int &x, int &y) { 69 if (!b) 70 d = a, x = 1, y = 0; 71 else { 72 exgcd(b, a % b, d, y, x); 73 y -= (a / b) * x; 74 } 75 } 76 77 int inv(int a, int n) { 78 int d, x, y; 79 exgcd(a, n, d, x, y); 80 return (x < 0) ? (x + n) : (x); 81 } 82 83 inline boolean init() { 84 return ~scanf("%d%d%d", &x, &m, &a) && (x || m || a); 85 } 86 87 int cs; 88 HashMap mp; 89 inline int ind(int pro, int x, int a, int p) { 90 mp.clear(); 91 cs = sqrt(p - 1 + 0.5); 92 int ainv = inv(x, p), iap = a * 1ll * qpow(ainv, cs - 1, p) % p; 93 for (int i = cs - 1; ~i; i--, iap = iap * 1ll * x % p) 94 mp.insert(iap, i); 95 int cp = qpow(x, cs, p), pw = pro; 96 for (int i = 0; i < p; i += cs, pw = pw * 1ll * cp % p) 97 if (~mp[pw]) 98 return mp[pw] + i; 99 return -1; 100 } 101 102 int exind(int x, int a, int m) { 103 if (a == 1) return 0; 104 int d, k = 0, pro = 1; 105 while ((d = gcd(x, m)) != 1) { 106 if (a % d) return -1; 107 if (pro == a) return k; 108 a /= d, m /= d, pro = (pro * 1ll * (x / d)) % m, k++; 109 } 110 int rt = ind(pro, x % m, a, m); 111 return (~rt) ? (rt + k) : (-1); 112 } 113 114 inline void solve() { 115 int res = exind(x % m, a % m, m); 116 if (res == -1) 117 puts("No Solution"); 118 else 119 printf("%d\n", res); 120 } 121 122 int main() { 123 while (init()) 124 solve(); 125 return 0; 126 }
