∑MathAdvanced

Discrete Logarithm

Key Points

•
The discrete logarithm problem asks for x such that $g^{x}$ ≡ h (mod p) in a multiplicative group modulo a prime p.
•
Exponentiation modulo p is fast, but inverting it (finding x) is generally hard and underpins cryptographic security.
•
Baby-step giant-step solves discrete logs in O( $p$ ) time and O( $p$ ) memory by meeting in the middle.
•
Pollard’s rho for discrete log also runs in expected O( $p$ ) time but uses only O(1) memory by searching for collisions in a pseudo-random walk.
•
Pohlig–Hellman is very fast when p−1 is smooth (factors into small primes) by reducing to small-subgroup discrete logs and combining with the CRT.
•
Solutions exist only if h lies in the subgroup generated by g, and the answer is unique modulo the order of g.
•
For competitive programming, BSGS handles primes up to about 10^12 comfortably; Pollard’s rho helps when memory is tight.
•
Care is needed with modular inverses modulo p−1, handling non-primitive g, and preventing 64-bit overflow in modular multiplication.

Prerequisites

→Modular arithmetic and congruences — Understanding equivalence modulo p, multiplicative inverses, and exponent reduction is essential for defining and solving DLP.
→Fast modular exponentiation — All algorithms require efficient computation of g^k mod p in O(log k) time.
→Hash maps and time–space trade-offs — BSGS relies on storing and looking up baby steps using a hash table.
→Extended Euclidean algorithm — Needed to compute modular inverses modulo p−1 and for CRT recombination.
→Chinese Remainder Theorem — Pohlig–Hellman recombines residues modulo prime powers into a single solution modulo p−1.
→Integer factorization basics — Pohlig–Hellman needs the prime-power factorization of p−1 to split the problem.
→Probability and the birthday paradox — Explains why Pollard’s rho finds collisions in O(√n) expected steps.
→C++ 64-bit arithmetic and overflow control — Implementations must avoid overflow in modular multiplication using 128-bit intermediates.

Detailed Explanation

Tap terms for definitions

01Overview

The discrete logarithm problem (DLP) is a fundamental task in number theory and cryptography: given a prime p, a base g, and a value h in the multiplicative group modulo p, find an integer x such that g^x ≡ h (mod p). In the cyclic group (Z/pZ)^\times of size n = p − 1, exponentiation is easy using fast powering, but finding the exponent from the result appears hard for general parameters. This asymmetry is what gives security to protocols like Diffie–Hellman key exchange, ElGamal encryption, and various signature schemes. From an algorithmic standpoint, the DLP invites clever strategies: meet-in-the-middle (baby-step giant-step) achieves O(\sqrt{n}) time and memory; Pollard’s rho reduces memory to O(1) with similar expected time using collision search; and Pohlig–Hellman exploits special structure when n factors into small primes. In programming contests, you’ll often meet DLP with manageable primes, where these algorithms are practical with 64-bit arithmetic. In cryptographic settings, parameters are chosen so that these algorithms are still infeasible, making DLP a hard problem in general.

02Intuition & Analogies

Think of g as a button that, when pressed once, moves you along a circular track of size p − 1 by a fixed but hidden stride; pressing it x times lands you at h. Exponentiation is like repeatedly pressing the button—easy. The discrete log is the reverse: standing at h, how many presses did it take to get here? Brute force would be trying 0, 1, 2, ... until you match h—too slow. Baby-step giant-step is like arranging a meeting: you take lots of small baby steps from the start (store where you land after j presses), and separately take giant leaps backwards from h (undo big chunks of presses) until you land somewhere you’ve seen. When the paths cross, you reconstruct the exact number of presses. Pollard’s rho imagines two people wandering the circle with rules depending on where they are; because the space is finite, their paths eventually cycle and collide (birthday paradox). From the collision, you deduce how many button presses separate them, revealing x. Finally, if the circle’s size breaks into small loops (p − 1 has small prime factors), Pohlig–Hellman cracks the problem piece by piece: solve how many presses modulo each small loop (small subgroup), then stitch those answers together with the Chinese Remainder Theorem to get the full count. These ideas turn a seemingly impossible reversal into tractable computations—provided the parameters are friendly.

03Formal Definition

Let p be prime and consider the multiplicative group G = (

Z

Z

)^{

\times

}, which is cyclic of order

n = p

− 1. Given g

\in

G and h

\in

G, the discrete logarithm problem asks for x

\in

Z

such that

g^{x}

\equiv

(mod p)

. If g has order ord(g) =

n_{g}

(a divisor of n), then a solution exists if and only if h

\in

⟨

⟩

, and the solution is unique modulo

n_{g}

. Algorithms: (1) Baby-step giant-step (Shanks) sets m =

⌈

n_{g}

⌉

, precomputes

g^{j}

for 0

\leq

j < m

, and searches for i such that h

g^{- im}

matches a baby step; time and memory are both O(

n_{g}

). (2) Pollard’s rho for DLP defines a pseudo-random iteration on states

X = g^{a}

h^{b}

and uses cycle-detection; a collision

X = X^{'}

yields

g^{a - a^{'}} = h^{b^{'} - b}

, from which x is solved modulo

n_{g}

g cd

(b'−b,

n_{g}

) = 1 (or handled more generally). Expected time is O(

n_{g}

), memory O(1). (3) Pohlig–Hellman factors

n_{g}

\prod

i^{e_{i}}

, solves x modulo each q_

i^{e_{i}}

by lifting discrete logs in subgroups of order

q_{i}

, then combines with the Chinese Remainder Theorem. Correctness relies on group theory (cyclic subgroups, orders) and number theory (modular inverses, CRT).

04When to Use

Use baby-step giant-step when you need a deterministic, relatively straightforward O(\sqrt{n}) algorithm and can afford O(\sqrt{n}) memory—typical for p up to around 10^{12} in competitive programming. Use Pollard’s rho for DLP when memory is tight or when a randomized algorithm is acceptable; it often runs faster in practice with much less memory and is well-suited for larger p where storing \sqrt{p} entries is impossible. Apply Pohlig–Hellman when p − 1 is smooth (i.e., factors into small primes); then the problem decomposes into several small discrete logs that are fast to solve and recombine, offering dramatic speedups. In mixed cases, factor p − 1 partially: if it has a large smooth part and a single large prime factor r, solve the smooth part with Pohlig–Hellman and the remaining modulo r by BSGS or Pollard’s rho. In cryptography, you typically do not solve DLPs—rather, systems choose p and g so that all known sub-exponential attacks (in prime fields, the best general algorithms are sub-exponential but not polynomial) are still infeasible at the chosen sizes. In contests, always consider constraints: pick BSGS for moderate sizes, switch to Pollard’s rho when memory is limited, and try Pohlig–Hellman if p − 1 factors nicely.

⚠️Common Mistakes

• Assuming g is a primitive root without checking: if ord(g) < p − 1 and h ∉ ⟨g⟩, no solution exists. Either ensure g is a generator or confirm solvability (e.g., factor p − 1 to get ord(g) and check membership). • Using Fermat inverses where not valid: inverses modulo p − 1 (for exponent arithmetic) require the extended Euclidean algorithm, not Fermat’s little theorem. • Overflow in modular multiplication: with 64-bit p, naive a*b % p overflows; use 128-bit intermediates (__int128) or Montgomery multiplication. • Incorrect BSGS step sizes: off-by-one in m = ⌈√n⌉ or building baby steps with duplicates can produce wrong answers. Store the smallest j for each value to avoid inflating answers. • Forgetting to reduce exponents modulo p − 1 when raising g^k (mod p): always compute k mod (p − 1) to keep numbers small and correct. • Pollard’s rho collision handling: failing to handle gcd(b2 − b1, n) > 1 leads to division by a non-invertible number; handle the general case or restart with new parameters. • Ignoring negative residues: when computing differences like (a1 − a2) mod n, normalize to [0, n − 1] before taking inverses. • Hash-map performance pitfalls: in BSGS, reserve() capacity to avoid rehash storms and consider custom hash for 64-bit if TLE. • Not validating the found x: always check g^x ≡ h (mod p) and, if needed, reduce x modulo ord(g) to canonical form.

Key Formulas

Discrete Logarithm Definition

g^{x} \equiv h (mod p)

Explanation: Find the exponent x such that raising g to x equals h modulo the prime p. The solution is unique modulo the order of g if it exists.

Group Order Mod Prime

n = ∣ (Z / p Z)^{\times} ∣ = p - 1

Explanation: The multiplicative group modulo a prime p has exactly p−1 nonzero elements, and it is cyclic. Exponent arithmetic is done modulo n.

BSGS Step Size

m = ⌈ n ⌉

Explanation: Choose m as the step size for baby and giant steps. This balances the two phases to achieve O(√n) time and space.

BSGS Matching Condition

g^{j} \equiv h (g^{- m})^{i} (mod p)

Explanation: Precompute all $g^{j}$ (baby steps) and search for i such that h·( $g^{- m}$ )^i hits one of them. Then $x = i \cdot m + j$ (mod n).

BSGS Complexity

T_{BSGS} (n) = O (n), S_{BSGS} (n) = O (n)

Explanation: Baby-step giant-step requires square-root time and memory relative to the group order n. This is feasible for moderate n.

Pollard’s Rho Collision Equation

g^{a} h^{b} = g^{a^{'}} h^{b^{'}} \Rightarrow g^{a - a^{'}} = h^{b^{'} - b}

Explanation: A collision between two states yields a linear congruence in x: (a−a') ≡ x (b'−b) (mod n). Solve it to recover x modulo n.

Solving the Collision

x \equiv (a - a^{'}) (b^{'} - b)^{- 1} (mod n)

Explanation: If gcd(b'−b, n)=1, invert (b'−b) modulo n and multiply to get x. If not, reduce by gcd and handle multiple solutions.

Factorization for Pohlig–Hellman

n = i = 1 \prod k q_{i}^{e_{i}}

Explanation: Decompose the group order into prime powers. The discrete log is solved modulo each q_ $i^{e_{i}}$ and recombined with CRT.

CRT Reconstruction

x \equiv x_{i} (mod q_{i}^{e_{i}}) \forall i \Rightarrow x \equiv i = 1 \sum k x_{i} M_{i} y_{i} (mod n)

Explanation: With $M_{i} = n$ / (q_ $i^{e_{i}}$ ) and $y_{i} = M_{i}^{- 1}$ mod q_ $i^{e_{i}}$ , this linear combination yields the unique solution modulo n.

Fermat’s Little Theorem

a^{p - 1} \equiv 1 (mod p)

Explanation: For prime p and a not divisible by p, exponent reduction modulo p−1 is valid: $a^{k}$ mod p depends only on k mod (p−1).

Pollard’s Rho Expected Time

E [steps of rho] = Θ (n)

Explanation: By the birthday paradox, a collision is expected after about √n random steps, giving the rho algorithm its time complexity.

Complexity Analysis

Let n denote the group order (for prime modulus p,

n = p

− 1). Baby-step giant-step (BSGS) sets m = ⌈√n⌉ and performs O(m) precomputations of

g^{j}

(baby steps) and up to O(m) giant steps of h·(

g^{- m}

)^i, with hash lookups per step. Hence time is O(√n) and memory is O(√n); the constant factors come from modular exponentiation and hash operations. Pollard’s rho for discrete log constructs a pseudo-random iteration over states

X = g^{a} h^{b}

and uses cycle detection. Under heuristic randomness, the expected number of steps to find a collision is Θ(√n) (birthday bound). Each step is O(1) modular multiplications, so the expected time is O(√n) with O(1) extra memory; however, restarts may be needed if the linear congruence is not solvable due to gcd issues, slightly increasing constants. Pohlig–Hellman’s complexity depends on the factorization n = ∏ q_

i^{e_{i}}

. It reduces the problem to k prime-power subproblems; solving modulo

q_{i}

typically costs O(√

q_{i}

) using BSGS (or O(

q_{i}

) by brute force if

q_{i}

is tiny), and lifting to q_

i^{e_{i}}

repeats this O(

e_{i}

) times. Overall, the cost is roughly O(∑

e_{i}

√{

q_{i}

}) modular exponentiations plus CRT combination in O(k lo

g^{2}

n). If n is B-smooth, this can be far below √n. Space for Pohlig–Hellman is dominated by the largest subgroup discrete log (O(√max

q_{i}

)) plus negligible CRT storage. In practice, algorithm choice trades memory for runtime and depends critically on n’s factorization and input sizes.

Code Examples

Baby-step Giant-step (BSGS) for discrete log modulo prime p

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 using int64 = long long;
5 using i128 = __int128_t;
6 
7 // Safe modular multiplication (avoids 64-bit overflow)
8 static inline int64 mod_mul(int64 a, int64 b, int64 mod) {
9     return (int64)((i128)a * b % mod);
10 }
11 
12 // Fast modular exponentiation: computes a^e mod mod
13 static int64 mod_pow(int64 a, long long e, int64 mod) {
14     int64 r = 1 % mod;
15     int64 x = (a % mod + mod) % mod;
16     while (e > 0) {
17         if (e & 1) r = mod_mul(r, x, mod);
18         x = mod_mul(x, x, mod);
19         e >>= 1;
20     }
21     return r;
22 }
23 
24 // Modular inverse modulo prime p using Fermat (p must be prime and a != 0)
25 static int64 mod_inv_prime(int64 a, int64 p) {
26     return mod_pow(a, p - 2, p);
27 }
28 
29 // Baby-step Giant-step: solve g^x ≡ h (mod p), p prime (so group order n = p-1)
30 // Returns x in [0, n-1] if exists, else -1
31 int64 discrete_log_bsgs(int64 p, int64 g, int64 h) {
32     if (p == 2) return h % p; // trivial
33     g %= p; h %= p;
34     if (h == 1 % p) return 0; // g^0 = 1
35     if (g == 0) return (h == 0 ? 1 : -1); // outside multiplicative group
36 
37     const int64 n = p - 1; // order of multiplicative group modulo prime p
38 
39     int64 m = (int64)ceil(sqrt((long double)n));
40 
41     // Baby steps: store g^j -> j for j in [0, m)
42     unordered_map<int64, int64> table;
43     table.reserve((size_t)m * 2);
44     table.max_load_factor(0.7);
45 
46     int64 baby = 1 % p;
47     for (int64 j = 0; j < m; ++j) {
48         // only store first occurrence to keep smallest j
49         if (!table.count(baby)) table[baby] = j;
50         baby = mod_mul(baby, g, p);
51     }
52 
53     // Compute factor = g^{-m} mod p
54     int64 g_inv = mod_inv_prime(g, p);
55     int64 factor = mod_pow(g_inv, m, p);
56 
57     // Giant steps: h * (g^{-m})^i
58     int64 giant = h % p;
59     for (int64 i = 0; i <= m; ++i) {
60         auto it = table.find(giant);
61         if (it != table.end()) {
62             // x = i*m + j (mod n)
63             int64 x = i * m + it->second;
64             x %= n;
65             // verify
66             if (mod_pow(g, x, p) == h) return x;
67         }
68         giant = mod_mul(giant, factor, p);
69     }
70     return -1; // no solution (likely h not in <g>)
71 }
72 
73 int main() {
74     ios::sync_with_stdio(false);
75     cin.tie(nullptr);
76 
77     // Example usage:
78     // Input: p g h  (p prime)
79     // Output: x or -1 if no solution
80     long long p, g, h;
81     if (!(cin >> p >> g >> h)) return 0;
82 
83     long long x = discrete_log_bsgs(p, g, h);
84     cout << x << "\n";
85     return 0;
86 }
87

This program implements Shanks’ baby-step giant-step algorithm to solve g^x ≡ h (mod p) for prime p. It precomputes baby steps g^j for j < m = ⌈√(p−1)⌉ and then walks over giant steps h·(g^{−m})^i, searching for a collision via a hash map. It uses 128-bit intermediates to avoid overflow and Fermat’s little theorem for modular inverses modulo p. If a collision is found, it reconstructs x and verifies the solution.

Time: O(\sqrt{p})Space: O(\sqrt{p})

Pollard’s rho for discrete logarithm (prime modulus, O(1) memory)

1 #include <bits/stdc++.h>
2 using namespace std;
3 using int64 = long long;
4 using i128 = __int128_t;
5 
6 static inline int64 mod_mul(int64 a, int64 b, int64 mod) { return (int64)((i128)a * b % mod); }
7 static int64 mod_pow(int64 a, uint64_t e, int64 mod) {
8     int64 r = 1 % mod, x = (a % mod + mod) % mod;
9     while (e) { if (e & 1) r = mod_mul(r, x, mod); x = mod_mul(x, x, mod); e >>= 1; }
10     return r;
11 }
12 
13 // Extended GCD: returns gcd(a,b) and finds x,y with ax+by=g
14 static long long exgcd(long long a, long long b, long long &x, long long &y) {
15     if (b == 0) { x = (a >= 0 ? 1 : -1); y = 0; return llabs(a); }
16     long long x1, y1; long long g = exgcd(b, a % b, x1, y1);
17     x = y1; y = x1 - (a / b) * y1; return g;
18 }
19 
20 static long long mod_inv(long long a, long long mod) {
21     a %= mod; if (a < 0) a += mod;
22     long long x, y; long long g = exgcd(a, mod, x, y);
23     if (g != 1) return -1; // not invertible
24     x %= mod; if (x < 0) x += mod; return x;
25 }
26 
27 struct State { long long X, a, b; };
28 
29 // Partitioned update function for Pollard's Rho (3-way split)
30 struct Stepper {
31     long long p, g, h, n; // n = p-1
32     Stepper(long long p_, long long g_, long long h_) : p(p_), g(g_%p_), h(h_%p_), n(p_-1) {}
33     inline void step(State &s) const {
34         long long cls = s.X % 3;
35         if (cls == 0) { // multiply by h: (a, b+1)
36             s.X = mod_mul(s.X, h, p);
37             s.b = (s.b + 1) % n;
38         } else if (cls == 1) { // square: (2a, 2b)
39             s.X = mod_mul(s.X, s.X, p);
40             s.a = (2 * s.a) % n;
41             s.b = (2 * s.b) % n;
42         } else { // multiply by g: (a+1, b)
43             s.X = mod_mul(s.X, g, p);
44             s.a = (s.a + 1) % n;
45         }
46     }
47 };
48 
49 // Pollard's rho for DLP: solve g^x = h (mod p), p prime.
50 // Returns x in [0, n-1] or -1 if failure.
51 long long pollards_rho_dlp(long long p, long long g, long long h, int max_restarts = 10) {
52     if (p == 2) return h % p;
53     g %= p; h %= p; if (h == 1 % p) return 0;
54     long long n = p - 1;
55 
56     Stepper stp(p, g, h);
57     std::mt19937_64 rng((uint64_t)chrono::high_resolution_clock::now().time_since_epoch().count());
58 
59     for (int attempt = 0; attempt < max_restarts; ++attempt) {
60         // Random seed state: X = g^a h^b
61         long long a0 = (long long)(rng() % n);
62         long long b0 = (long long)(rng() % n);
63         State tort{ mod_mul(mod_pow(g, a0, p), mod_pow(h, b0, p), p), a0, b0 };
64         State hare = tort;
65 
66         // Floyd's cycle detection
67         for (uint64_t iter = 1; iter <= (uint64_t)5 * (uint64_t)sqrt((long double)n) + 1000000ULL; ++iter) {
68             stp.step(tort);
69             stp.step(hare); stp.step(hare);
70             if (tort.X == hare.X) {
71                 // g^{a1} h^{b1} = g^{a2} h^{b2}
72                 long long a1 = tort.a, b1 = tort.b;
73                 long long a2 = hare.a, b2 = hare.b;
74                 long long r = (a1 - a2) % n; if (r < 0) r += n;
75                 long long s = (b2 - b1) % n; if (s < 0) s += n;
76                 if (s == 0) break; // failure; restart
77                 long long d = std::gcd(s, n);
78                 if (r % d != 0) break; // inconsistent; restart
79                 long long n2 = n / d;
80                 long long s2 = s / d;
81                 long long r2 = r / d;
82                 long long inv = mod_inv(s2 % n2, n2);
83                 if (inv == -1) break; // not invertible; restart
84                 long long x0 = (i128)r2 * inv % n2;
85                 // There are d solutions: x ≡ x0 (mod n2). Try the smallest.
86                 long long x = x0;
87                 if (mod_pow(g, x, p) == h) return x;
88                 // Try other coset solutions if needed
89                 for (long long k = 1; k < d && k < 5; ++k) {
90                     long long cand = (x0 + k * n2) % n;
91                     if (mod_pow(g, cand, p) == h) return cand;
92                 }
93                 break; // give up this attempt
94             }
95         }
96     }
97     return -1;
98 }
99 
100 int main() {
101     ios::sync_with_stdio(false);
102     cin.tie(nullptr);
103 
104     long long p, g, h; // p prime
105     if (!(cin >> p >> g >> h)) return 0;
106     long long x = pollards_rho_dlp(p, g, h);
107     cout << x << "\n";
108     return 0;
109 }
110

This program implements Pollard’s rho algorithm for discrete log in a prime field. It performs a pseudo-random walk over states X = g^a h^b, partitions the state space into three classes, and uses Floyd’s cycle detection to find a collision. The collision yields a linear congruence in x modulo n = p − 1, which the code solves using the extended Euclidean algorithm (to handle non-invertible cases). It restarts with random seeds on failure and verifies candidate solutions.

Time: O(\sqrt{p}) expectedSpace: O(1)

Pohlig–Hellman (with subgroup BSGS) for smooth p−1 and CRT recombination

1 #include <bits/stdc++.h>
2 using namespace std;
3 using int64 = long long;
4 using i128 = __int128_t;
5 
6 static inline int64 mod_mul(int64 a, int64 b, int64 mod) { return (int64)((i128)a * b % mod); }
7 static int64 mod_pow(int64 a, unsigned long long e, int64 mod) {
8     int64 r = 1 % mod, x = (a % mod + mod) % mod;
9     while (e) { if (e & 1ULL) r = mod_mul(r, x, mod); x = mod_mul(x, x, mod); e >>= 1ULL; }
10     return r;
11 }
12 
13 // Extended GCD for inverses modulo composite moduli (e.g., p-1 or CRT steps)
14 static long long exgcd(long long a, long long b, long long &x, long long &y) {
15     if (b == 0) { x = (a >= 0 ? 1 : -1); y = 0; return llabs(a); }
16     long long x1, y1; long long g = exgcd(b, a % b, x1, y1);
17     x = y1; y = x1 - (a / b) * y1; return g;
18 }
19 static long long mod_inv_ll(long long a, long long mod) {
20     a %= mod; if (a < 0) a += mod; long long x, y; long long g = exgcd(a, mod, x, y);
21     if (g != 1) return -1; x %= mod; if (x < 0) x += mod; return x;
22 }
23 
24 // BSGS restricted to subgroup of known order 'ord': solve base^x = target (mod p), 0<=x<ord
25 static long long bsgs_known_order(long long p, long long base, long long target, long long ord) {
26     base %= p; target %= p; if (target == 1 % p) return 0;
27     long long m = (long long)ceil(sqrt((long double)ord));
28     unordered_map<long long, long long> table; table.reserve((size_t)m * 2); table.max_load_factor(0.7);
29     long long baby = 1 % p;
30     for (long long j = 0; j < m; ++j) {
31         if (!table.count(baby)) table[baby] = j;
32         baby = mod_mul(baby, base, p);
33     }
34     long long inv_base = mod_inv_ll(base % p, p); if (inv_base == -1) return -1;
35     long long factor = 1 % p;
36     // factor = base^{-m}
37     for (long long i = 0; i < m; ++i) factor = mod_mul(factor, inv_base, p);
38 
39     long long giant = target;
40     for (long long i = 0; i <= m; ++i) {
41         auto it = table.find(giant);
42         if (it != table.end()) {
43             long long x = i * m + it->second;
44             if (x < ord && mod_pow(base, x, p) == target) return x;
45         }
46         giant = mod_mul(giant, factor, p);
47     }
48     return -1;
49 }
50 
51 // Trial division factorization of n into prime powers (sufficient for many contest limits)
52 static vector<pair<long long,int>> factorize_pp(long long n) {
53     vector<pair<long long,int>> fac;
54     for (long long p = 2; p * p <= n; ++p) {
55         if (n % p == 0) {
56             int e = 0; while (n % p == 0) { n /= p; ++e; }
57             fac.push_back({p, e});
58         }
59     }
60     if (n > 1) fac.push_back({n, 1});
61     return fac;
62 }
63 
64 // Pohlig–Hellman for prime modulus p: solve g^x = h (mod p), assuming n = p-1 is reasonably smooth.
65 // Returns x in [0, n-1] or -1 if fails.
66 long long pohlig_hellman(long long p, long long g, long long h) {
67     g %= p; h %= p; if (h == 1 % p) return 0;
68     long long n = p - 1; // group order
69 
70     // Factor n into prime powers
71     auto fac = factorize_pp(n);
72 
73     long long x = 0; // will hold solution modulo n
74     long long M = 1; // current modulus
75 
76     for (auto [q, e] : fac) {
77         // Solve x modulo q^e
78         long long qe = 1; for (int i = 0; i < e; ++i) qe *= q;
79         long long E0 = n / qe;
80         long long g0 = mod_pow(g, E0, p); // order q^e
81         long long h0 = mod_pow(h, E0, p);
82 
83         long long xi = 0; // digits in base q
84         long long qpow = 1; // q^j
85         for (int j = 0; j < e; ++j) {
86             // E = n / q^{j+1}
87             long long Ej = n;
88             for (int t = 0; t < j + 1; ++t) Ej /= q; // careful but fast enough for small e
89 
90             // c_j = h^{E} * (g^{xi})^{−E} mod p
91             // Reduce exponent modulo p-1 to keep it small
92             long long gx = mod_pow(g, (unsigned long long)(xi % (p - 1)), p);
93             long long c1 = mod_pow(h, (unsigned long long)Ej, p);
94             long long inv = mod_pow(gx, (unsigned long long)(p - 2), p);
95             long long c = mod_mul(c1, inv, p);
96 
97             // base for this digit: a = g^{E} has order q
98             long long a = mod_pow(g, (unsigned long long)Ej, p);
99             long long dj = bsgs_known_order(p, a, c, q);
100             if (dj == -1) return -1; // failure
101             xi += dj * qpow;
102             qpow *= q;
103         }
104         // Now x ≡ xi (mod q^e). Merge with CRT into current (x mod M).
105         // Since q^e are pairwise coprime across different primes, we can use standard CRT.
106         long long m1 = M, a1 = x;
107         long long m2 = qe, a2 = xi % qe;
108         // Solve: x ≡ a1 (mod m1), x ≡ a2 (mod m2)
109         long long inv = mod_inv_ll(m1 % m2, m2);
110         if (inv == -1) return -1; // should not happen as m1 and m2 are coprime
111         long long t = (a2 - a1) % m2; if (t < 0) t += m2;
112         long long k = (i128)t * inv % m2;
113         x = a1 + k * m1;
114         M = m1 * m2;
115         x %= M;
116     }
117 
118     // Verify
119     if (mod_pow(g, (unsigned long long)(x % (p - 1)), p) == h) return x % (p - 1);
120     return -1;
121 }
122 
123 int main() {
124     ios::sync_with_stdio(false);
125     cin.tie(nullptr);
126 
127     long long p, g, h; // p prime
128     if (!(cin >> p >> g >> h)) return 0;
129     long long x = pohlig_hellman(p, g, h);
130     cout << x << "\n";
131     return 0;
132 }
133

This program implements Pohlig–Hellman for prime modulus p, which exploits the factorization of n = p − 1. It decomposes the problem into prime-power subgroups, solves each digit using a small-subgroup discrete log (via BSGS with known order q), lifts through the exponent e, and combines the results using the Chinese Remainder Theorem. It uses trial division for factoring (sufficient for many contest settings) and verifies the final answer.

Time: O(\sum e_i \sqrt{q_i}) for subgroup logs + CRT, where p−1 = \prod q_i^{e_i}Space: O(\sqrt{\max q_i})

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	using int64 = long long;
5	using i128 = __int128_t;
6
7	// Safe modular multiplication (avoids 64-bit overflow)
8	static inline int64 mod_mul(int64 a, int64 b, int64 mod) {
9	return (int64)((i128)a * b % mod);
10	}
11
12	// Fast modular exponentiation: computes a^e mod mod
13	static int64 mod_pow(int64 a, long long e, int64 mod) {
14	int64 r = 1 % mod;
15	int64 x = (a % mod + mod) % mod;
16	while (e > 0) {
17	if (e & 1) r = mod_mul(r, x, mod);
18	x = mod_mul(x, x, mod);
19	e >>= 1;
20	}
21	return r;
22	}
23
24	// Modular inverse modulo prime p using Fermat (p must be prime and a != 0)
25	static int64 mod_inv_prime(int64 a, int64 p) {
26	return mod_pow(a, p - 2, p);
27	}
28
29	// Baby-step Giant-step: solve g^x ≡ h (mod p), p prime (so group order n = p-1)
30	// Returns x in [0, n-1] if exists, else -1
31	int64 discrete_log_bsgs(int64 p, int64 g, int64 h) {
32	if (p == 2) return h % p; // trivial
33	g %= p; h %= p;
34	if (h == 1 % p) return 0; // g^0 = 1
35	if (g == 0) return (h == 0 ? 1 : -1); // outside multiplicative group
36
37	const int64 n = p - 1; // order of multiplicative group modulo prime p
38
39	int64 m = (int64)ceil(sqrt((long double)n));
40
41	// Baby steps: store g^j -> j for j in [0, m)
42	unordered_map<int64, int64> table;
43	table.reserve((size_t)m * 2);
44	table.max_load_factor(0.7);
45
46	int64 baby = 1 % p;
47	for (int64 j = 0; j < m; ++j) {
48	// only store first occurrence to keep smallest j
49	if (!table.count(baby)) table[baby] = j;
50	baby = mod_mul(baby, g, p);
51	}
52
53	// Compute factor = g^{-m} mod p
54	int64 g_inv = mod_inv_prime(g, p);
55	int64 factor = mod_pow(g_inv, m, p);
56
57	// Giant steps: h * (g^{-m})^i
58	int64 giant = h % p;
59	for (int64 i = 0; i <= m; ++i) {
60	auto it = table.find(giant);
61	if (it != table.end()) {
62	// x = i*m + j (mod n)
63	int64 x = i * m + it->second;
64	x %= n;
65	// verify
66	if (mod_pow(g, x, p) == h) return x;
67	}
68	giant = mod_mul(giant, factor, p);
69	}
70	return -1; // no solution (likely h not in <g>)
71	}
72
73	int main() {
74	ios::sync_with_stdio(false);
75	cin.tie(nullptr);
76
77	// Example usage:
78	// Input: p g h (p prime)
79	// Output: x or -1 if no solution
80	long long p, g, h;
81	if (!(cin >> p >> g >> h)) return 0;
82
83	long long x = discrete_log_bsgs(p, g, h);
84	cout << x << "\n";
85	return 0;
86	}
87