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

View this file on GitHub
Last update: 2026-06-20 09:18:49+09:00
Include: #include "math/prime_factorization.hpp"

Overview

Deterministic Miller-Rabin primality testing and Pollard-Rho factorization for the full uint64_t range. This is the general-purpose choice when values are too large for a sieve.

For example,

\[360 = 2^3 \cdot 3^2 \cdot 5.\]

prime_factors(360) returns 2, 2, 2, 3, 3, 5, while prime_factorize(360) returns the pairs (2, 3), (3, 2), and (5, 1).

Use PrimeSieve instead when every query is bounded by a reasonably small known limit. A sieve has a setup cost but answers many small queries faster.

Euler’s Totient Function

Euler’s totient function $\varphi(n)$ counts the integers from 1 through n that are coprime to n. Two integers are coprime when their greatest common divisor is 1.

For example, the integers coprime to 12 are

1, 5, 7, 11

so $\varphi(12) = 4$.

If the distinct prime divisors of n are $p_1, p_2, \ldots, p_k$, then

\[\varphi(n) = n \left(1-\frac{1}{p_1}\right) \left(1-\frac{1}{p_2}\right) \cdots \left(1-\frac{1}{p_k}\right).\]

The products in this displayed formula are multiplicative: each parenthesized factor is multiplied with n and with the other factors.

Equivalently, and often more clearly for implementation:

\[\varphi(n) = n \prod_{p \mid n}\frac{p-1}{p}.\]

Common uses include:

counting reduced fractions with denominator n;
working with multiplicative groups modulo n;
reducing exponents using Euler’s theorem: $a^{\varphi(n)} \equiv 1 \pmod n$ when a and n are coprime.

The library names this function euler_phi in this header and totient in PrimeSieve.

Mobius Function

The Mobius function $\mu(n)$ is defined from the prime factorization of n:

$\mu(1) = 1$;
$\mu(n) = 0$ if some prime square divides n;
otherwise, $\mu(n) = (-1)^k$, where k is the number of distinct prime factors.

Examples:

Value	Factorization	Mobius value	Reason
`1`	empty product	`1`	Definition
`6`	`2 * 3`	`1`	Two distinct prime factors
`30`	`2 * 3 * 5`	`-1`	Three distinct prime factors
`12`	`2^2 * 3`	`0`	A prime square divides it

Its main competitive-programming use is inclusion-exclusion over divisors. Mobius inversion says that if

\[F(n) = \sum_{d \mid n} f(d),\]

then

\[f(n) = \sum_{d \mid n} \mu(d) F(n/d).\]

This often converts counts over divisors into counts with an exact gcd, or counts all pairs into counts of coprime pairs.

The conventional spelling is “Möbius”; the API uses ASCII name mobius.

How the Algorithms Fit Together

Miller-Rabin tests whether a number is prime without trying every possible divisor. For 64-bit integers, the fixed witness set used here makes the result deterministic rather than merely probable.

Pollard-Rho searches for a nontrivial divisor of a composite number using a pseudo-random sequence and gcd computations. Once it finds a divisor, the implementation recursively factors both pieces and uses Miller-Rabin to know when a piece is already prime.

API

bool is_prime(uint64_t value);

std::vector<uint64_t> prime_factors(uint64_t value);

std::vector<std::pair<uint64_t, int>> prime_factorize(uint64_t value);

std::vector<uint64_t> divisors(uint64_t value);

uint64_t euler_phi(uint64_t value);

int mobius(uint64_t value);

Inputs use uint64_t, so negative integers are not accepted. prime_factorize stores each prime as uint64_t and its exponent as int. The Mobius function returns int because its result is always -1, 0, or 1; the other numeric result uses uint64_t.

Function	Description
`is_prime(x)`	Deterministically tests whether `x` is prime.
`prime_factors(x)`	Returns prime factors with multiplicity in increasing order.
`prime_factorize(x)`	Returns `(prime, exponent)` pairs in increasing order.
`divisors(x)`	Returns all positive divisors in increasing order.
`euler_phi(x)`	Returns Euler’s totient function.
`mobius(x)`	Returns the Mobius function.

All functions except is_prime require x >= 1.

divisors(x) includes both 1 and x. For example, the divisors of 12 are 1, 2, 3, 4, 6, 12.

Complexity

Miller-Rabin uses a fixed seven-base witness set and takes $O(\log x)$ modular multiplications.

Pollard-Rho has probabilistic expected running time of roughly $O(x^{1/4})$ for finding a factor, and is very fast for ordinary 64-bit competitive-programming inputs. The returned result is deterministic in content even though the search uses pseudo-random polynomial parameters.

Example

#include "math/prime_factorization.hpp"

#include <cstdint>
#include <iostream>

int main() {
    uint64_t value = 360;
    for (const auto& factor : m1une::math::prime_factorize(value)) {
        std::cout << factor.first << "^" << factor.second << "\n";
    }

    std::cout << m1une::math::euler_phi(12) << "\n";  // 4
    std::cout << m1une::math::mobius(30) << "\n";     // -1
}

Required by

Verified with

Code

#ifndef M1UNE_MATH_PRIME_FACTORIZATION_HPP
#define M1UNE_MATH_PRIME_FACTORIZATION_HPP 1

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

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

#endif  // M1UNE_MATH_PRIME_FACTORIZATION_HPP

#line 1 "math/prime_factorization.hpp"



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

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

m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

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

Overview

Euler’s Totient Function

Mobius Function

How the Algorithms Fit Together

API

Complexity

Example

Required by

Verified with

Code

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

Overview

Euler’s Totient Function

Mobius Function

How the Algorithms Fit Together

API

Complexity

Example

Required by

Verified with

Code

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