Primitive Root
(math/primitive_root.hpp)

• View this file on GitHub
• Last update: 2026-06-24 15:25:00+09:00
• Include: #include "math/primitive_root.hpp"

Overview

This header finds a primitive root modulo mod.

A primitive root is a residue whose powers generate every invertible residue modulo mod. Equivalently, its multiplicative order is Euler’s totient function $\varphi(\text{mod})$.

Primitive roots do not exist for every modulus. They exist exactly for

• 2;
• 4;
• p^k, where p is an odd prime and k >= 1;
• 2 * p^k, where p is an odd prime and k >= 1.

API

bool has_primitive_root(uint64_t mod);

uint64_t primitive_root(uint64_t mod);

has_primitive_root(mod) returns whether a primitive root exists.

primitive_root(mod) requires mod >= 2. It returns the smallest positive primitive root modulo mod, or 0 when no primitive root exists.

The return value 0 is only used as a sentinel: a primitive root is always coprime to mod and lies in [1, mod).

Complexity

The modulus and $\varphi(\text{mod})$ are factored with Pollard-Rho, so this supports the full uint64_t range.

After factorization, candidates are tested from small to large. Each candidate uses one modular exponentiation for each distinct prime divisor of $\varphi(\text{mod})$.

Example

#include "math/primitive_root.hpp"

#include <cstdint>
#include <iostream>

int main() {
    uint64_t mod = 998244353;
    uint64_t root = m1une::math::primitive_root(mod);
    std::cout << root << "\n";  // 3

    std::cout << m1une::math::primitive_root(8) << "\n";  // 0
}

Depends on

• 64-bit Prime Factorization (math/prime_factorization.hpp)

Required by

• Math All (math/all.hpp)
• Multidimensional Convolution (math/multivariate_convolution.hpp)

Verified with

• verify/math/math_algorithms.test.cpp
• verify/math/multivariate_convolution_cyclic.test.cpp
• verify/math/multivariate_convolution_truncated.test.cpp
• verify/math/primitive_root.test.cpp

Code

#ifndef M1UNE_MATH_PRIMITIVE_ROOT_HPP
#define M1UNE_MATH_PRIMITIVE_ROOT_HPP 1

#include <cassert>
#include <cstdint>
#include <numeric>
#include <utility>
#include <vector>

#include "prime_factorization.hpp"

namespace m1une {
namespace math {

inline bool has_primitive_root(uint64_t mod) {
    if (mod == 2 || mod == 4) return true;
    if (mod < 2) return false;

    uint64_t odd_part = mod;
    if ((odd_part & 1) == 0) {
        odd_part >>= 1;
        if ((odd_part & 1) == 0) return false;
    }

    return prime_factorize(odd_part).size() == 1;
}

// Returns the smallest positive primitive root modulo mod.
// Returns 0 when no primitive root exists.
inline uint64_t primitive_root(uint64_t mod) {
    assert(mod >= 2);
    if (mod == 2) return 1;
    if (!has_primitive_root(mod)) return 0;

    const uint64_t phi = euler_phi(mod);
    const std::vector<std::pair<uint64_t, int>> factors = prime_factorize(phi);
    for (uint64_t candidate = 2; candidate < mod; candidate++) {
        if (std::gcd(candidate, mod) != 1) continue;

        bool generator = true;
        for (const auto& factor : factors) {
            if (internal::power_mod(candidate, phi / factor.first, mod) == 1) {
                generator = false;
                break;
            }
        }
        if (generator) return candidate;
    }
    return 0;
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_PRIMITIVE_ROOT_HPP

#line 1 "math/primitive_root.hpp"



#include <cassert>
#include <cstdint>
#include <numeric>
#include <utility>
#include <vector>

#line 1 "math/prime_factorization.hpp"



#include <algorithm>
#line 10 "math/prime_factorization.hpp"

namespace m1une {
namespace math {

namespace internal {

inline uint64_t multiply_mod(uint64_t a, uint64_t b, uint64_t mod) {
    return static_cast<uint64_t>(static_cast<unsigned __int128>(a) * b % mod);
}

inline uint64_t power_mod(uint64_t base, uint64_t exponent, uint64_t mod) {
    uint64_t result = 1;
    while (exponent > 0) {
        if (exponent & 1) result = multiply_mod(result, base, mod);
        base = multiply_mod(base, base, mod);
        exponent >>= 1;
    }
    return result;
}

inline uint64_t pollard_random() {
    static uint64_t state = 0x123456789abcdef0ULL;
    state += 0x9e3779b97f4a7c15ULL;
    uint64_t value = state;
    value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL;
    value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL;
    return value ^ (value >> 31);
}

}  // namespace internal

inline bool is_prime(uint64_t value) {
    if (value < 2) return false;
    for (uint64_t prime : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (value % prime == 0) return value == prime;
    }

    uint64_t odd_part = value - 1;
    int power_of_two = 0;
    while ((odd_part & 1) == 0) {
        odd_part >>= 1;
        power_of_two++;
    }

    for (uint64_t base : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
        if (base % value == 0) continue;
        uint64_t x = internal::power_mod(base % value, odd_part, value);
        if (x == 1 || x == value - 1) continue;

        bool composite = true;
        for (int i = 1; i < power_of_two; i++) {
            x = internal::multiply_mod(x, x, value);
            if (x == value - 1) {
                composite = false;
                break;
            }
        }
        if (composite) return false;
    }
    return true;
}

namespace internal {

inline uint64_t pollard_rho(uint64_t value) {
    for (uint64_t prime : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (value % prime == 0) return prime;
    }

    while (true) {
        const uint64_t constant = pollard_random() % (value - 1) + 1;
        uint64_t y = pollard_random() % (value - 1) + 1;
        uint64_t x = 0;
        uint64_t saved_y = 0;
        uint64_t gcd = 1;
        uint64_t segment_length = 1;

        auto advance = [&](uint64_t current) {
            return static_cast<uint64_t>(
                (static_cast<unsigned __int128>(multiply_mod(current, current, value)) + constant) % value);
        };

        while (gcd == 1) {
            x = y;
            for (uint64_t i = 0; i < segment_length; i++) y = advance(y);

            for (uint64_t offset = 0; offset < segment_length && gcd == 1; offset += 128) {
                saved_y = y;
                uint64_t product = 1;
                const uint64_t block = std::min<uint64_t>(128, segment_length - offset);
                for (uint64_t i = 0; i < block; i++) {
                    y = advance(y);
                    const uint64_t difference = x > y ? x - y : y - x;
                    product = multiply_mod(product, difference, value);
                }
                gcd = std::gcd(product, value);
            }
            segment_length <<= 1;
        }

        if (gcd == value) {
            do {
                saved_y = advance(saved_y);
                const uint64_t difference = x > saved_y ? x - saved_y : saved_y - x;
                gcd = std::gcd(difference, value);
            } while (gcd == 1);
        }
        if (gcd != value) return gcd;
    }
}

inline void factor_recursively(uint64_t value, std::vector<uint64_t>& factors) {
    if (value == 1) return;
    if (is_prime(value)) {
        factors.push_back(value);
        return;
    }
    const uint64_t divisor = pollard_rho(value);
    factor_recursively(divisor, factors);
    factor_recursively(value / divisor, factors);
}

}  // namespace internal

inline std::vector<uint64_t> prime_factors(uint64_t value) {
    assert(value >= 1);
    std::vector<uint64_t> result;
    internal::factor_recursively(value, result);
    std::sort(result.begin(), result.end());
    return result;
}

inline std::vector<std::pair<uint64_t, int>> prime_factorize(uint64_t value) {
    std::vector<uint64_t> factors = prime_factors(value);
    std::vector<std::pair<uint64_t, int>> result;
    for (uint64_t prime : factors) {
        if (result.empty() || result.back().first != prime) {
            result.emplace_back(prime, 1);
        } else {
            result.back().second++;
        }
    }
    return result;
}

inline std::vector<uint64_t> divisors(uint64_t value) {
    std::vector<uint64_t> result = {1};
    for (const auto& factor : prime_factorize(value)) {
        const int current_size = int(result.size());
        uint64_t power = 1;
        for (int exponent = 1; exponent <= factor.second; exponent++) {
            power *= factor.first;
            for (int i = 0; i < current_size; i++) {
                result.push_back(result[i] * power);
            }
        }
    }
    std::sort(result.begin(), result.end());
    return result;
}

inline uint64_t euler_phi(uint64_t value) {
    assert(value >= 1);
    uint64_t result = value;
    for (const auto& factor : prime_factorize(value)) {
        result = result / factor.first * (factor.first - 1);
    }
    return result;
}

inline int mobius(uint64_t value) {
    assert(value >= 1);
    int result = 1;
    for (const auto& factor : prime_factorize(value)) {
        if (factor.second >= 2) return 0;
        result = -result;
    }
    return result;
}

}  // namespace math
}  // namespace m1une


#line 11 "math/primitive_root.hpp"

namespace m1une {
namespace math {

inline bool has_primitive_root(uint64_t mod) {
    if (mod == 2 || mod == 4) return true;
    if (mod < 2) return false;

    uint64_t odd_part = mod;
    if ((odd_part & 1) == 0) {
        odd_part >>= 1;
        if ((odd_part & 1) == 0) return false;
    }

    return prime_factorize(odd_part).size() == 1;
}

// Returns the smallest positive primitive root modulo mod.
// Returns 0 when no primitive root exists.
inline uint64_t primitive_root(uint64_t mod) {
    assert(mod >= 2);
    if (mod == 2) return 1;
    if (!has_primitive_root(mod)) return 0;

    const uint64_t phi = euler_phi(mod);
    const std::vector<std::pair<uint64_t, int>> factors = prime_factorize(phi);
    for (uint64_t candidate = 2; candidate < mod; candidate++) {
        if (std::gcd(candidate, mod) != 1) continue;

        bool generator = true;
        for (const auto& factor : factors) {
            if (internal::power_mod(candidate, phi / factor.first, mod) == 1) {
                generator = false;
                break;
            }
        }
        if (generator) return candidate;
    }
    return 0;
}

}  // namespace math
}  // namespace m1une

Back to top page