ix. Primitive Roots and Discrete Log

The order ${\text{ord}}_n(a)$ is the least $k>0$ with $a^k\equiv 1(\mathrm{mod}n)$ ($\gcd (a,n)=1$); it divides $\phi (n)$ and the Carmichael exponent $\lambda (n)$. A primitive root has order $\phi (n)$; it exists iff $n\in \{1,2,4,p^k,2p^k\}$ ($p$ odd prime), and then ${\mathbb{Z}}_n^{\ast }$ is cyclic: every unit is $g^x$ and multiplication becomes addition of exponents (indices). Discrete log ($a^x\equiv b$) and discrete roots ($x^k\equiv a$) both reduce to it.

Task	Method	Cost
order of $a$	factor $\phi (n)$ (or $\lambda (n)$), strip primes while $a^{k/q}\equiv 1$	$O({\log }^{2}n)$ after factoring
primitive root mod $p$	test $g=2,3,\dots$: $g^{(p-1)/q}\ne 1$ for all primes $q\mid p-1$	small $g$ suffice in practice ($O({\log }^{6}p)$ under GRH)
number of primitive roots	$\phi (\phi (n))$	—
$a^x\equiv b(\mathrm{mod}m)$, any $m$	baby-step giant-step (+ gcd reduction for non-coprime $a$)	$O(\sqrt{m})$
same, $p-1$ smooth	Pohlig–Hellman	$O(\sum e_i(\log p+\sqrt{q_i}))$
$x$ known to lie in $[l,r]$	Pollard kangaroo	$O(\sqrt{r-l})$, $O(1)$ memory
$x^k\equiv a(\mathrm{mod}p)$	$x=g^y$, solve $ky\equiv \text{ind}(a)(\mathrm{mod}p-1)$	BSGS + linear congruence

Order, primitive roots, discrete log, roots

vector<ll> primeFactorsOf(ll n) { vector<ll> f; for (ll d = 2; d * d <= n; d++) if (n % d == 0) { f.push_back(d); while (n % d == 0) n /= d; } if (n > 1) f.push_back(n); return f; }
ll primitiveRootPrime(ll p) {                              // p prime
    if (p == 2) return 1;
    auto qs = primeFactorsOf(p - 1);
    for (ll g = 2;; g++) { bool ok = true; for (ll q : qs) if (power(g, (p - 1) / q, p) == 1) { ok = false; break; } if (ok) return g; }
}
ll discreteLog(ll a, ll b, ll m) {                         // smallest x >= 0 with a^x = b (mod m), -1 if none; any m >= 1
    a %= m; b %= m; if (m == 1) return 0;
    ll k = 1 % m, add = 0, g;
    while ((g = std::gcd(a, m)) > 1) {                     // peel common factors: a^x = b  ->  (a/g) a^{x-1} = b/g (mod m/g)
        if (b == k) return add;
        if (b % g) return -1;
        b /= g, m /= g, add++, k = (i128)k * (a / g) % m;
    }
    ll n = (ll)sqrtl(m) + 1, an = 1;                       // now solve k * a^y = b with gcd(a, m) = 1
    unordered_map<ll, ll> vals; vals.reserve(2 * n);
    for (ll q = 0, cur = b; q <= n; q++) { vals[cur] = q; cur = (i128)cur * a % m; }     // baby steps b a^q (larger q overwrites)
    for (ll i = 0; i < n; i++) an = (i128)an * a % m;
    for (ll p = 1, cur = k; p <= n; p++) {
        cur = (i128)cur * an % m;
        if (auto it = vals.find(cur); it != vals.end()) return n * p - it->second + add;
    }
    return -1;
}
vector<ll> discreteRoots(ll k, ll a, ll p) {               // all x in [0, p) with x^k = a (mod p), p prime
    a %= p;
    if (k == 0) { vector<ll> v; if (a == 1 % p) { v.resize(p); iota(all(v), 0); } return v; }   // x^0 = 1 (0^0 = 1)
    if (a == 0) return {0};
    ll g = primitiveRootPrime(p), ia = discreteLog(g, a, p), n = p - 1;
    ll d = std::gcd(k, n); if (ia % d) return {};
    ll step = n / d, y0 = (ll)((i128)(ia / d) * invMod((k / d) % step, step) % step);   // k y = ia (mod n)
    vector<ll> res; for (ll t = 0; t < d; t++) res.push_back(power(g, y0 + t * step, p));
    sort(all(res)); return res;
}

The general-purpose library below adds multiplicativeOrder, carmichaelLambda, jacobi, pohligHellman (fast when $p-1$ is smooth) and pollardPm1 (finds a factor $d$ of $n$ when some prime $q\mid n$ has smooth $q-1$), with its own 64-bit factorisation.

// ---------- basic number theory (number-theory/basic-number-theory.cpp) ----------
u64 mulMod(u64 a, u64 b, u64 m) {
    return (u64)((u128)a * b % m);
}
 
u64 powMod(u64 a, u64 e, u64 m) {
    u64 res = 1 % m;
    a %= m;
    for (; e > 0; e >>= 1) {
        if (e & 1) res = mulMod(res, a, m);
        a = mulMod(a, a, m);
    }
    return res;
}
 
bool isPrimeTrial(i64 n) {
    if (n < 4) return n >= 2;
    if (n % 2 == 0 || n % 3 == 0) return false;
    for (i64 i = 5; i * i <= n; i += 6) {
        if (n % i == 0 || n % (i + 2) == 0) return false;
    }
    return true;
}
 
// prime factorization as (prime, exponent) pairs in increasing order of primes.
vector<pair<i64, int>> factorizeTrial(i64 n) {
    vector<pair<i64, int>> res;
    for (i64 p = 2; p * p <= n; p += (p == 2 ? 1 : 2)) {
        if (n % p != 0) continue;
        int e = 0;
        while (n % p == 0) n /= p, e++;
        res.push_back({p, e});
    }
    if (n > 1) res.push_back({n, 1});
    return res;
}
 
// all positive divisors of n >= 1 in increasing order.
vector<i64> divisorsTrial(i64 n) {
    vector<i64> small, large;
    for (i64 i = 1; i * i <= n; i++) {
        if (n % i != 0) continue;
        small.push_back(i);
        if (i != n / i) large.push_back(n / i);
    }
    small.insert(small.end(), large.rbegin(), large.rend());
    return small;
}
 
i64 eulerPhi(i64 n) {
    i64 res = n;
    for (auto [p, e] : factorizeTrial(n)) res -= res / p;
    return res;
}
 
template <class T>
vector<T> divisorsFromFactors(const vector<pair<T, int>>& factors) {
    vector<T> res = {1};
    for (auto [p, e] : factors) {
        int cnt = res.size();
        T pw = 1;
        for (int i = 0; i < e; i++) {
            pw *= p;
            for (int j = 0; j < cnt; j++) res.push_back(res[j] * pw);
        }
    }
    sort(res.begin(), res.end());
    return res;
}
 
// ---------- Miller-Rabin and Pollard-Rho factorization (number-theory/factorization.cpp) ----------
bool isPrime(u64 n) {
    static constexpr u64 bases[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    if (n < 2) return false;
    for (u64 p : bases) {
        if (n % p == 0) return n == p;
    }
    int s = __builtin_ctzll(n - 1);
    u64 d = (n - 1) >> s;
    for (u64 a : bases) {
        u64 x = powMod(a, d, n);
        if (x == 1 || x == n - 1) continue;
        int i = 1;
        for (; i < s; i++) {
            x = mulMod(x, x, n);
            if (x == n - 1) break;
        }
        if (i == s) return false;
    }
    return true;
}
 
u64 pollardRho(u64 n) {
    if (n % 2 == 0) return 2;
    static mt19937_64 gen(chrono::steady_clock::now().time_since_epoch().count());
    const u64 batch = 128;
    u64 c = 1;
    auto step = [&](u64 x) {   // x^2 + c mod n without overflow
        u64 r = mulMod(x, x, n) + c;
        return (r < c || r >= n) ? r - n : r;
    };
    auto dist = [](u64 a, u64 b) { return a > b ? a - b : b - a; };
    while (true) {
        c = gen() % (n - 1) + 1;
        u64 y = gen() % n, x = y, ys = y, q = 1, g = 1;
        for (u64 len = 1; g == 1; len <<= 1) {
            x = y;
            for (u64 i = 0; i < len; i++) y = step(y);
            for (u64 k = 0; k < len && g == 1; k += batch) {
                ys = y;
                for (u64 i = 0; i < batch && i < len - k; i++) {
                    y = step(y);
                    q = mulMod(q, dist(x, y), n);
                }
                g = gcd(q, n);
            }
        }
        if (g == n) {   // the batch overshot: redo it one step at a time
            do {
                ys = step(ys);
                g = gcd(dist(x, ys), n);
            } while (g == 1);
        }
        if (g != n) return g;
    }
}
 
vector<pair<u64, int>> factorize(u64 n) {
    assert(n >= 1);
    vector<u64> primes;
    for (u64 p = 2; p < 100 && p * p <= n; p++) {
        while (n % p == 0) primes.push_back(p), n /= p;
    }
    auto rec = [&](auto&& self, u64 m) -> void {
        if (m == 1) return;
        if (isPrime(m)) {
            primes.push_back(m);
            return;
        }
        u64 d = pollardRho(m);
        self(self, d);
        self(self, m / d);
    };
    rec(rec, n);
    sort(primes.begin(), primes.end());
    vector<pair<u64, int>> res;
    for (u64 p : primes) {
        if (!res.empty() && res.back().first == p) res.back().second++;
        else res.push_back({p, 1});
    }
    return res;
}
 
vector<u64> divisors(u64 n) {
    return divisorsFromFactors(factorize(n));
}
 
// ---------- multiplicative order, Carmichael function, Jacobi symbol, Pohlig-Hellman, Pollard p - 1 ----------
u64 multiplicativeOrder(u64 a, u64 n) {
    if (n == 1) return 1;
    u64 phi = n;
    for (auto [p, e] : factorize(n)) phi = phi / p * (p - 1);
    u64 ord = phi;
    for (auto [q, e] : factorize(phi)) {
        while (ord % q == 0 && powMod(a, ord / q, n) == 1) ord /= q;
    }
    return ord;
}
 
u64 carmichaelLambda(u64 n) {
    u64 res = 1;
    for (auto [p, e] : factorize(n)) {
        u64 pk = 1;
        for (int i = 0; i < e - 1; i++) pk *= p;
        u64 l = pk * (p - 1);   // phi(p^e)
        if (p == 2 && e >= 3) l /= 2;
        res = res / gcd(res, l) * l;
    }
    return res;
}
 
int jacobi(i64 a, i64 n) {
    a %= n;
    if (a < 0) a += n;
    int res = 1;
    while (a != 0) {
        while (a % 2 == 0) {
            a /= 2;
            if (n % 8 == 3 || n % 8 == 5) res = -res;
        }
        swap(a, n);
        if (a % 4 == 3 && n % 4 == 3) res = -res;
        a %= n;
    }
    return n == 1 ? res : 0;
}
 
i64 pohligHellman(u64 g, u64 h, u64 p) {
    g %= p, h %= p;
    if (p == 2) return h == 1 ? 0 : (g == 1 ? -1 : 1);
    u64 order = multiplicativeOrder(g, p);
    if (powMod(h, order, p) != 1) return -1;   // h not in <g>
    u64 x = 0, m = 1;   // x mod m, combined by CRT
    for (auto [q, e] : factorize(order)) {
        u64 qe = 1;
        for (int i = 0; i < e; i++) qe *= q;
        u64 gq = powMod(g, order / q, p);   // element of order q
        // baby-step giant-step table for logarithms base gq (group of size q)
        u64 step = (u64)sqrtl((long double)q) + 1;
        unordered_map<u64, u64> table;
        u64 cur = 1;
        for (u64 j = 0; j < step; j++) table.emplace(cur, j), cur = mulMod(cur, gq, p);
        u64 giant = powMod(gq, q - step % q, p);   // gq^(-step)
        auto logQ = [&](u64 y) -> i64 {
            u64 v = y;
            for (u64 i = 0; i <= step; i++) {
                auto it = table.find(v);
                if (it != table.end()) return (i64)((i * step + it->second) % q);
                v = mulMod(v, giant, p);
            }
            return -1;
        };
        // digits of x mod q^e
        u64 xk = 0, qk = 1;
        u64 ginv = powMod(g, order - 1, p);
        for (int k = 0; k < e; k++) {
            u64 hk = mulMod(h, powMod(ginv, xk, p), p);
            hk = powMod(hk, order / (qk * q), p);
            i64 d = logQ(hk);
            if (d < 0) return -1;
            xk += d * qk, qk *= q;
        }
        // combine x = xk mod qe with the previous congruence
        auto inverse = [](u64 a, u64 mod) {   // a^-1 mod `mod`, gcd(a, mod) = 1
            i64 t = 0, newT = 1, r = mod, newR = a % mod;
            while (newR) {
                i64 q2 = r / newR;
                tie(t, newT) = pair(newT, t - q2 * newT);
                tie(r, newR) = pair(newR, r - q2 * newR);
            }
            return (u64)((t % (i64)mod + (i64)mod) % (i64)mod);
        };
        u64 t = (u64)((i128)((xk % qe + qe - x % qe) % qe) * inverse(m % qe, qe) % qe);
        x += m * t, m *= qe;
    }
    return x % order;
}
 
u64 pollardPm1(u64 n, u64 bound = 100000) {
    if (n % 2 == 0) return n == 2 ? 0 : 2;
    u64 a = 2;
    for (u64 k = 2; k <= bound; k++) {
        a = powMod(a, k, n);
        if (k % 64 == 0 || k == bound) {
            u64 d = gcd(a == 0 ? n : a - 1, n);
            if (d == n) return 0;
            if (d > 1) return d;
        }
    }
    return 0;
}

Proof: order divides, primitive-root test, existence and count

Order divides: if $a^k\equiv 1$ write $k=q\text{ord}+r$; then $a^r\equiv 1$ with $r<\text{ord}$ forces $r=0$. So $\text{ord}\mid \phi (n)$ (Euler) and $\text{ord}$ can be found by starting from $\phi (n)$ and dividing out primes while the power stays 1. Test: $g$ has order $p-1$ iff no proper divisor of $p-1$ is a period, iff $g^{(p-1)/q}\ne 1$ for each prime $q\mid p-1$ (every proper divisor divides some $(p-1)/q$). Existence mod $p$: a polynomial of degree $d$ has at most $d$ roots in ${\mathbb{F}}_p$, so for each $d\mid p-1$ there are at most $\phi (d)$ elements of order $d$ (they are roots of $x^d-1$ generated by any one of them); since ${\sum }_{d\mid p-1}\phi (d)=p-1$ every bound is tight, in particular $\phi (p-1)\ge 1$ elements have order $p-1$. Lifting to $p^k$ and $2p^k$ and the failure for other $n$ (a product of two groups of even order is not cyclic) are standard; in a cyclic group of order $m$ the generators are $g^j$ with $\gcd (j,m)=1$, giving $\phi (m)$ of them.

Proof: baby-step giant-step, the gcd reduction, and discrete roots

BSGS: for coprime $a$, write $x=np-q$ with $1\le p\le n$, $0\le q\le n$ (every $x\in [0,n^2]$ has such a form, and $n^2\ge m>$ the order, so some solution is in range if any exists). Then $k{a}^{np}\equiv b{a}^q$, which is found by looking up giant steps among stored baby steps; scanning $p$ upward and keeping the largest $q$ for each value returns the smallest $x$. Reduction: if $g=\gcd (a,m)>1$ and $x\ge 1$, then $a^x\equiv b(\mathrm{mod}m)$ requires $g\mid b$ and is equivalent to $\frac{a}{g}{a}^{x-1}\equiv \frac{b}{g}(\mathrm{mod}\frac{m}{g})$; the factor $\frac{a}{g}$ is accumulated in $k$ and each step divides $m$ by at least 2, so after $O(\log m)$ steps $a$ is coprime to $m$; the check b == k catches solutions $x$ smaller than the number of peeled steps. Roots: with $x=g^y$ and $a=g^{\text{ind}(a)}$, $x^k\equiv a⟺ky\equiv \text{ind}(a)(\mathrm{mod}p-1)$, a linear congruence with $\gcd (k,p-1)$ solutions when solvable.

Proof sketch: Pohlig–Hellman and Pollard's kangaroo

Pohlig–Hellman: if $p-1=\prod q_i^{e_i}$, then $x\mathrm{mod}{q}_i^{e_i}$ is recovered digit by digit in base $q_i$: raising $g^x=h$ to the power $(p-1)/q_i^{j+1}$ isolates the $j$-th digit as a discrete log in the subgroup of order $q_i$ (solved by BSGS in $O(\sqrt{q_i})$); CRT combines the residues. Kangaroo: a tame kangaroo jumps from $g^r$ with pseudo-random jumps depending on the current value and leaves a trap; a wild one starts from $h$ and follows the same deterministic rule, so once it lands on a point of the tame path it follows it into the trap. With mean jump $\approx \sqrt{r-l}/2$, the expected work is $O(\sqrt{r-l})$ with constant memory (heuristic, birthday-style analysis).