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

扩展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 }
ex-BSGS

 

posted @ 2018-02-27 13:54  阿波罗2003  阅读(...)  评论(... 编辑 收藏