Baby Step Giant Step 算法

2015.02.26 12:47 Thu| 20 visits oi_2015| 2015_算法笔记| Text

Baby Step Giant Step 算法

这个算法真是起了一个好名字。这么浪漫竟然被安到了一个数论算法头上,真是……

可以解出如下方程 $A^x \equiv B \mod P$ 的最小非负整数解,其中 P 是素数。

这个算法的原理实在是非常简单:

令答案 $x = i * m + j$,其中 $m = \sqrt P,\quad i, j \le m$。于是得到 $A^x = A^{i * m} * A ^ j$。这样就可以分块解决这个问题了:

先计算出对于 $0 \le i \le m$ 的 $A^i \mod P$ 并插入到哈希表(map)中。之后枚举 $i$,解方程 $x *A^{i * m}=B$。若解得的 x 在哈希表中出现过,就得到了答案。

下面贴上我的不知为何常数巨大的代码(想要常数?去玩手写 hash 咯!):

#include <bits/stdc++.h>
using namespace std;

int power(int a, int b, int p)
{
    int re = 1;
    while (b)
    {
        if (b & 1) re = 1ll * re * a % p;
        a = 1ll * a * a % p; b >>= 1;
    }
    return re;
}

int bsgs(int a, int b, int p)
{
    map<int, int> m;
    a %= p;
    if (!a) return b ? -1 : 1;
    int block = ceil(sqrt(p));
    int t = 1; m[t] = 0;
    for (int i = 1; i < block; ++i)
    {
        t = 1ll * t * a % p;
        if (!m[t]) m[t] = i;
    }
    int inv = power(a, p - block - 1, p);
    for (int i = 0; i < block; ++i)
    {
        if (m.find(b) != m.end())
            return i * block + m[b];
        b = 1ll * b * inv % p;
    }
    return -1;
}

int main()
{
    long long a, b, p, ans;
    scanf("%d%d%d", &a, &b, &p);
    ans = bsgs(a, b, p);
    if (ans == -1)
        puts("no solution");
    else printf("%d\n", ans);;
    return 0;
}

附赠果题一双

BZOJ3239 Discrete Logging

BZOJ2242 [SDOI2011]计算器

分为三个子问题,第一个是简简单单快速幂,第二个是简简单单 EXGCD,第三个是简简单单 BSGS。

Code

#include <bits/stdc++.h>
using namespace std;

inline int read()
{
    int x = 0, f = 1; char ch = getchar();
    while (ch < '0' || ch > '9') { if (ch == '-') f = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') { x = x * 10 + ch - '0'; ch = getchar(); }
    return x * f;
}

long long power(long long a, long long b, long long p)
{
    long long re = 1;
    while (b)
    {
        if (b & 1) re = re * a % p;
        a = a * a % p; b >>= 1;
    }
    return re;
}

pair<long long, long long> exgcd(long long a, long long b)
{
    if (b == 0) return make_pair(1, 0);
    pair<long long, long long> t = exgcd(b, a % b);
    return make_pair(t.second, t.first - a / b * t.second);
}

long long bsgs(long long a, long long b, long long p)
{
    map<long long, int> m;
    a %= p;
    if (!a) return b ? -1 : 1;
    int block = ceil(sqrt(p));
    long long t = 1; m[t] = 0;
    for (int i = 1; i < block; ++i)
    {
        t = t * a % p;
        if (!m[t]) m[t] = i;
    }
    long long inv = power(a, p - block - 1, p);
    for (int i = 0; i < block; ++i)
    {
        if (m.find(b) != m.end())
            return i * block + m[b];
        b = b * inv % p;
    }
    return -1;
}

int main()
{
    int t = read(), k = read();
    while (t--)
    {
        long long y = read(), z = read(), p = read(), ans;
        if (k == 1) cout << power(y, z, p) << '\n';
        else
        {
            if (k == 2)
                ans = (exgcd(y, p).first % p * z % p + p) % p,
                ans = ans ? ans : -1;
            else ans = bsgs(y, z, p);
            if (ans == -1)
                cout << "Orz, I cannot find x!\n";
            else cout << ans << '\n';
        }
    }
    return 0;
}

Extended Baby Step Giant Step 算法

传统的 BSGS 算法,很不幸,不能用来解决当 P 不是素数时的情况。

考虑 G 同时是 A 和 C 的因数。若 G 不是 B 的因数,可以确定无解。否则令 $B' = B / G, C' = C / G$。假设 x 不为 0,有 $(A/G') * A^{x-1} \equiv B' \mod C'$。

重复令这个约数为 $\gcd(A, C)$,经过几次这样的操作后,A 和 C 将会互素。设这种操作一共进行 cnt 次,每一次的约数之积为 T,则 $A^{cnt}/T * A^{x-cnt} \equiv B' mod C'$。到这里就可以套用 BSGS 算法解出 x 了。

似乎并没有结束,不要忘记上述解法解出的 x 不小于 cnt,对于判断 x 是否小于 cnt,解决办法就是枚举 $\log P$ 之内的答案暴力检验。

代码(写了好久,开始时乱七八糟)

int gcd(int a, int b)
{
    if (b == 0) return a;
    return gcd(b, a % b);
}

pair<int, int> exgcd(int a, int b)
{
    if (b == 0) return make_pair(1, 0);
    pair<int, int> t = exgcd(b, a % b);
    return make_pair(t.second, t.first - a / b * t.second);
}

int exbsgs(int a, int b, int p)
{
    if (p == 1) return 0;
    int t = 1, d = 1, cnt = 0;
    a %= p;
    for (int i = 0; i < 32; ++i)
    {
        if (t == b) return i;
        t = 1ll * t * a % p;
    }
    while ((t = gcd(a, p)) != 1)
    {
        if (b % t) return -1;
        b /= t; p /= t; d = 1ll * d * a / t % p;
        ++cnt;
    }
    map<int, int> m;
    int block = ceil(sqrt((double)p));
    m[t = 1] = 0;
    for (int i = 1; i < block; ++i)
    {
        t = 1ll * t * a % p;
        if (!m[t]) m[t] = i;
    }
    int k = power(a, block, p);
    for (int i = 0; i < block; ++i)
    {
        t = (1ll * exgcd(d, p).first * b % p + p) % p;
        if (m.find(t) != m.end())
            return i * block + m[t] + cnt;
        d = 1ll * d * k % p;
    }
    return -1;
}