Baby Step Giant Step 及其扩展算法

对于求解 A^x ≡ B (mod C) (0 <= x < C)这样的问题，可以使用 Baby Step Giant Step及其扩展算法解决。

对于BSGS，AC大神在他的博客中有一篇很详细的介绍，相信大部分人都是看他的文章学会BSGS的。这里我再做出一点整理。

Baby Step Giant Step

原始的BSGS只能解决C为素数的情况。

设m = ceil(sqrt(C)), x = i * m + j，那么A^x = (A^m)^i * A^j, (0 <= i < m, 0 <= j < m)。然后可以枚举i，O(sqrt(C))级别的枚举。

对于一个枚举出来的i，令D = (A^m)^i现在问题转化为求D * A^j ≡ B (mod C)。如果把A^j当作一个整体，那么套上exgcd就可以解出来了（而且因为C是质数，A是C的倍数的情况容易特判，除此之外必有(D, C) = 1，所以一定有解）

求出了A^j，现在的问题就是我怎么知道j是多少？AC大神给出了一个非常聪明的方法，先用O(sqrt(C))的时间，将A^j全部存进hash表里面。然后只要查表就在O(1)的时间内知道j是多少了。

Extended Baby Step Giant Step

扩展过的BSGS不要求C为素数。大致的做法和没扩展的一样，只是有些修改。

因为同余有一条性质

令 d = gcd(A, C), A = ad, B = bd, C = cd
则     ad ≡ bd (mod cd)
等价于 a  ≡ b  (mod c )

因此，开始前，我们先执行消因子：

D = 1, cnt = 0
while gcd(A, C) != 1:
    if (B % gcd(A, C) != 0)
        No Solution
    B /= gcd(A, C)
    C /= gcd(A, C)
    D *= A / gcd(A, C)
    cnt += 1

执行完了以后，问题变成了求D * A^(x-cnt) ≡ B (mod C)，令x = i * m + j + cnt，后面的做法就跟BSGS一样了。

但是注意到i, j >= 0，因此x >= cnt，但是我们明显存在<= cnt的情况，所以在消因子之前要做一次O(log(C))的枚举，直接验证A^i % C == B?，这样就避免了这种情况。

题目

BSGS：

  • poj2417

Extended BSGS：

  • poj3243
  • hdu2815
  • zju3254
  • spoj3105 MOD

BSGS Code Example

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

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>

typedef long long Int64;
class Hash
{
  static const int HASHMOD = 7679977;
  int top, hash[HASHMOD+100], value[HASHMOD+100], stack[1<<16];
  int locate(const int x) const
  {
    int h = x % HASHMOD;
    while (hash[h] != -1 && hash[h] != x) ++h;
    return h;
  }
public:
  Hash(): top(0) { memset(hash, 0xFF, sizeof(hash)); }
  void insert(const int x, const int v)
  {
    const int h = locate(x);
    if (hash[h] == -1)
      hash[h] = x, value[h] = v, stack[++top] = h;
  }
  int get(const int x) const
  {
    const int h = locate(x);
    return hash[h] == x ? value[h] : -1;
  }
  void clear() { while (top) hash[stack[top--]] = -1; }
} hash;
struct Triple
{
  Int64 x, y, z;
  Triple(const Int64 a, const Int64 b, const Int64 c): x(a), y(b), z(c) { }
};
Triple ExtendedGCD(const Int64 a, const Int64 b)
{
  if (b == 0) return Triple(1, 0, a);
  const Triple last = ExtendedGCD(b, a%b);
  return Triple(last.y, last.x - a / b * last.y, last.z);
}
int A, B, C;
Int64 BabyStep(Int64 A, Int64 B, Int64 C)
{
  const int sqrtn = static_cast<int>(std::ceil(std::sqrt(C)));
  Int64 base = 1;
  hash.clear();
  for (int i = 0; i < sqrtn; ++i)
  {
    hash.insert(base, i);
    base = base*A % C;
  }
  Int64 i = 0, j = -1, D = 1;
  for (; i < sqrtn; ++i)
  {
    Triple res = ExtendedGCD(D, C);
    const int c = C / res.z;
    res.x = (res.x * B / res.z % c + c) % c;
    j = hash.get(res.x);
    if (j != -1) return i * sqrtn + j;
    D = D * base % C;
  }
  return -1;
}
int main()
{
  while (scanf("%d%d%d", &C, &A, &B) != EOF)
  {
    const Int64 ans = BabyStep(A, B, C);
    if (ans == -1)
      printf("no solution\n");
    else
      printf("%I64d\n", ans);
  }
}

Extended BSGS Code Example

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

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>

typedef long long Int64;
class Hash
{
  static const int HASHMOD = 3894229;
  int top, hash[HASHMOD+100], value[HASHMOD+100], stack[1<<16];
  int locate(const int x) const
  {
    int h = x % HASHMOD;
    while (hash[h] != -1 && hash[h] != x) ++h;
    return h;
  }
public:
  Hash(): top(0) { memset(hash, 0xFF, sizeof(hash)); }
  void insert(const int x, const int v)
  {
    const int h = locate(x);
    if (hash[h] == -1)
      hash[h] = x, value[h] = v, stack[++top] = h;
  }
  int get(const int x) const
  {
    const int h = locate(x);
    return hash[h] == x ? value[h] : -1;
  }
  void clear() { while (top) hash[stack[top--]] = -1; }
} hash;
struct Triple
{
  Int64 x, y, z;
  Triple() { }
  Triple(const Int64 a, const Int64 b, const Int64 c): x(a), y(b), z(c) { }
};
Triple ExtendedGCD(const int a, const int b)
{
  if (b == 0) return Triple(1, 0, a);
  const Triple last = ExtendedGCD(b, a%b);
  return Triple(last.y, last.x - a / b * last.y, last.z);
}
int ExtendedBabyStep(int A, int B, int C)
{
  Int64 tmp = 1, cnt = 0, D = 1;
  for (int i = 0; i < 32; ++i)
  {
    if (tmp == B) return i;
    tmp = tmp * A % C;
  }
  for (Triple res; (res = ExtendedGCD(A, C)).z != 1; ++cnt)
  {
    if (B % res.z) return -1;
    B /= res.z, C /= res.z;
    D = D * A/res.z % C;
  }
  const int sqrtn = static_cast<int>(std::ceil(std::sqrt(C)));
  hash.clear();
  Int64 base = 1;
  for (int i = 0; i < sqrtn; ++i)
  {
    hash.insert(base, i);
    base = base * A % C;
  }
  Int64 j = -1, i;
  for (i = 0; i < sqrtn; ++i)
  {
    Triple res = ExtendedGCD(D, C);
    const int c = C / res.z;
    res.x = (res.x * B/res.z % c + c) % c;
    j = hash.get(res.x);
    if (j != -1) return i * sqrtn + j + cnt;
    D = D * base % C;
  }
  return -1;
}

int main()
{
  for (int a, b, c; scanf("%d%d%d", &a, &c, &b) != EOF && (a || b || c); )
  {
    b %= c;
    const int ans = ExtendedBabyStep(a, b, c);
    if (ans == -1)
      printf("No Solution\n");
    else
      printf("%d\n", ans);
  }
}

Comments