BSGS Algorithm

The Baby-step Giant-step (BSGS) algorithm solves discrete logarithm problems of the form \(a^x \equiv n \pmod{p}\) for the smallest non-negative integer \(x\), where \(a\) and \(p\) are coprime constants.

Algorithm Overview

Express \(x\) as \(x = A \cdot \lceil \sqrt{p} \rceil - B\) where \(0 \leq A, B \leq \lceil \sqrt{p} \rceil\). This transforms the equation to:

\[a^{A \cdot \lceil \sqrt{p} \rceil} \equiv n \cdot a^B \pmod{p}\]

Precompute all possible \(n \cdot a^B\) values and store them in a hash table with their corresponding \(B\) values. Then iterate over possible \(A\) values to find a matching entry. The solution is given by \(x = A \cdot \lceil \sqrt{p} \rceil - B\). Time complexity is \(O(\sqrt{p})\), or \(O(\sqrt{p} \log p)\) with a std::map.

Implementation

#include <iostream>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;

ll power_mod(ll base, ll exp, ll mod) {
    ll result = 1;
    while (exp) {
        if (exp & 1) result = (result * base) % mod;
        base = (base * base) % mod;
        exp >>= 1;
    }
    return result;
}

ll bsgs(ll base, ll target, ll prime) {
    base %= prime, target %= prime;
    if (target == 1) return 0;
    if (base == 0) return (target == 0 ? 1 : -1);

    map<ll, ll> index_map;
    ll step_size = ceil(sqrt(prime));
    ll current = target;
    
    for (ll j = 0; j <= step_size; j++) {
        index_map[current] = j;
        current = (current * base) % prime;
    }

    ll base_step = power_mod(base, step_size, prime);
    current = base_step;
    for (ll i = 1; i <= step_size; i++) {
        if (index_map.count(current) && (i * step_size - index_map[current] >= 0))
            return i * step_size - index_map[current];
        current = (current * base_step) % prime;
    }
    return -1;
}

int main() {
    ll p, a, n;
    cin >> p >> a >> n;
    ll result = bsgs(a, n, p);
    if (result != -1) cout << result << endl;
    else cout << "no solution" << endl;
    return 0;
}

Extended BSGS

Extended BSGS (exBSGS) solves \(a^x \equiv n \pmod{p}\) when \(a\) and \(p\) are not necessarily coprime.

Algorithm Overview

Repeatedly divide the equation by \(d = \gcd(a, p)\) until \(a\) and \(p\) become coprime. After \(k\) iterations, the equation transforms to:

\[a^{x-k} \cdot \frac{a^k}{\prod d} \equiv \frac{n}{\prod d} \pmod{\frac{p}{\prod d}}\]

Solve the transformed equation using BSGS and add \(k\) to the solution.

Implementation

#include <iostream>
#include <map>
using namespace std;
typedef long long ll;

ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; }

ll bsgs(ll base, ll target, ll prime, ll coeff);

ll exbsgs(ll base, ll target, ll prime) {
    base %= prime, target %= prime;
    if (target == 1 || prime == 1) return 0;
    if (base == 0) return (target == 0 ? 1 : -1);

    ll iter_count = 0, accum = 1;
    while (true) {
        ll d = gcd(base, prime);
        if (d == 1) break;
        if (target % d != 0) return -1;
        iter_count++;
        prime /= d, target /= d;
        accum = (accum * (base / d)) % prime;
        if (accum == target) return iter_count;
    }
    ll result = bsgs(base, target, prime, accum);
    return (result == -1) ? -1 : result + iter_count;
}

int main() {
    ll a, p, n;
    while (cin >> a >> p >> n && a) {
        ll result = exbsgs(a, n, p);
        if (result != -1) cout << result << endl;
        else cout << "No Solution" << endl;
    }
    return 0;
}

Matrix BSGS

Matrix BSGS solves \(A^x \equiv B \pmod{p}\) for matrices \(A\) and \(B\), using matrix multiplication and hashing.

Implementation

#include <iostream>
#include <cmath>
#include <unordered_map>
using namespace std;

typedef unsigned long long ull;
const int MATRIX_SIZE = 70;
const ull HASH_BASE = 998244353;

struct Matrix {
    int size;
    int data[MATRIX_SIZE][MATRIX_SIZE];
    
    Matrix(int n = 0) : size(n) { 
        for (int i = 0; i < MATRIX_SIZE; i++)
            fill(data[i], data[i] + MATRIX_SIZE, 0);
    }

    ull hash_value() {
        ull val = 0;
        for (int i = 0; i < size; i++)
            for (int j = 0; j < size; j++)
                val = val * HASH_BASE + data[i][j];
        return val;
    }

    Matrix operator*(const Matrix& other) {
        Matrix res(size);
        for (int i = 0; i < size; i++)
            for (int k = 0; k < size; k++)
                for (int j = 0; j < size; j++)
                    res.data[i][j] = (res.data[i][j] + data[i][k] * other.data[k][j]) % mod;
        return res;
    }
};

int matrix_bsgs(Matrix base, Matrix target, int mod) {
    if (target.is_identity()) return 0;
    int step = sqrt(mod) + 1;
    Matrix current = target;
    unordered_map<ull, int> hash_map;
    Matrix base_step(base.size);
    base_step.set_identity();

    for (int j = 0; j < step; j++) {
        hash_map[current.hash_value()] = j;
        current = base * current;
        base_step = base_step * base;
    }

    Matrix temp = base_step;
    for (int i = 1; i <= step; i++) {
        ull hash_val = temp.hash_value();
        if (hash_map.count(hash_val)) return i * step - hash_map[hash_val];
        temp = base_step * temp;
    }
    return -1;
}