Number Theory
(math/number_theory.hpp)

View this file on GitHub
Last update: 2026-07-01 22:14:13+09:00
Include: #include "math/number_theory.hpp"

Overview

Fast integer number-theory primitives for the extended Euclidean algorithm, modular arithmetic, the Chinese remainder theorem, and floor sums. All functions support 64-bit inputs and use wider intermediates where products could overflow.

Extended Greatest Common Divisor

extended_gcd(a, b) returns (g, x, y) satisfying Bézout’s identity:

\[g = \gcd(a,b) = ax + by.\]

The returned gcd is nonnegative, and both inputs may be negative or zero. For a = b = 0, the result is (0, 0, 0). The coefficients are one valid solution; they are not otherwise normalized.

The nonnegative gcd must fit in long long. In particular, an input whose gcd is $2^{63}$ is outside the supported result range.

Congruences and Modular Arithmetic

The notation

\[x \equiv r \pmod m\]

means that x and r have the same remainder when divided by m. Equivalently, x - r is divisible by m.

For example, $23 \equiv 2 \pmod 7$ because both numbers leave remainder 2, and 23 - 2 = 21 is divisible by 7.

Modular Power

pow_mod(x, exponent, mod) computes

\[x^{\text{exponent}} \bmod \text{mod}\]

without constructing the enormous integer power first. Binary exponentiation uses only logarithmically many multiplications, so even exponents near $10^{18}$ are practical.

This is useful for modular counting answers, repeated transformations, and number-theory algorithms such as primality testing.

Modular Inverse

A modular inverse of x modulo mod is a number y satisfying

\[x y \equiv 1 \pmod{\text{mod}}.\]

For example, the inverse of 3 modulo 11 is 4, because $3 \cdot 4 = 12 \equiv 1 \pmod {11}$.

An inverse exists exactly when x and mod are coprime, meaning their greatest common divisor is 1. inv_mod(x, mod) requires this condition.

inv_gcd(x, mod) is the more general primitive. It returns (g, y), where g = gcd(x, mod) and

\[x y \equiv g \pmod{\text{mod}}.\]

For inv_gcd(6, 15), one valid result is (3, 3) because $6 \cdot 3 = 18 \equiv 3 \pmod {15}$.

Chinese Remainder Theorem

The Chinese remainder theorem combines several remainder constraints into one. For example:

\[\begin{aligned} x &\equiv 2 \pmod 3,\\ x &\equiv 3 \pmod 5,\\ x &\equiv 2 \pmod 7. \end{aligned}\]

The smallest nonnegative solution is 23. Every other solution is $23 + 105k$, because 105 is the least common multiple of the moduli. Therefore crt returns (23, 105).

The moduli do not need to be pairwise coprime. A solution merely needs to be consistent. For example, x cannot be both even and congruent to 1 modulo 4, so constraints x = 0 mod 2 and x = 1 mod 4 cause crt to return (0, 0).

Typical uses include:

combining answers computed under several moduli;
finding the first time that several periodic events coincide;
solving simultaneous divisibility or remainder constraints.

Floor Sum

floor_sum(n, mod, a, b) computes

\[\sum_{i=0}^{n-1} \left\lfloor \frac{a i+b}{\text{mod}} \right\rfloor.\]

For n = 4, mod = 10, a = 6, and b = 3, the terms are

\[\left\lfloor\frac{3}{10}\right\rfloor + \left\lfloor\frac{9}{10}\right\rfloor + \left\lfloor\frac{15}{10}\right\rfloor + \left\lfloor\frac{21}{10}\right\rfloor = 0 + 0 + 1 + 2 = 3.\]

A direct loop costs $O(n)$. This function repeatedly transforms the problem in a Euclidean-algorithm-like way and costs only logarithmic time.

Floor sums appear when:

counting lattice points below a line;
summing quotient or remainder patterns;
counting integers satisfying linear inequalities;
optimizing periodic arithmetic sequences.

API

std::tuple<long long, long long, long long> extended_gcd(
    long long a,
    long long b
);

long long pow_mod(
    long long x,
    unsigned long long exponent,
    long long mod
);

std::pair<long long, long long> inv_gcd(
    long long x,
    long long mod
);

long long inv_mod(long long x, long long mod);

std::pair<long long, long long> crt(
    const std::vector<long long>& remainders,
    const std::vector<long long>& moduli
);

long long floor_sum(
    long long n,
    long long mod,
    long long a,
    long long b
);

All scalar inputs and outputs use long long, except that pow_mod accepts an unsigned long long exponent. extended_gcd returns a three-element tuple; the other multi-value functions return std::pair<long long, long long>. crt reads its two vectors by const reference and does not modify them.

Function	Description	Complexity
`extended_gcd(a, b)`	Returns `(g, x, y)` with `g = gcd(a, b)` and `a * x + b * y = g`.	$O(\log(\max(	a	,	b	)+1))$
`pow_mod(x, exponent, mod)`	Returns `x` raised to `exponent` modulo `mod`.	$O(\log(\text{exponent}))$
`inv_gcd(x, mod)`	Returns the gcd and a normalized inverse of `x / gcd` modulo `mod / gcd`.	$O(\log(\text{mod}))$
`inv_mod(x, mod)`	Returns the modular inverse of `x`. The arguments must be coprime.	$O(\log(\text{mod}))$
`crt(remainders, moduli)`	Solves a possibly non-coprime system with the Chinese remainder theorem.	$O(k \log M)$
`floor_sum(n, mod, a, b)`	Returns the sum of `floor((a * i + b) / mod)` over `0 <= i < n`.	$O(\log(\max(a, mod)))$

crt returns the smallest nonnegative solution and the least common multiple of the moduli. It returns (0, 0) when no solution exists. The least common multiple must fit in long long.

floor_sum accepts negative a and b; its result must fit in long long.

Example

#include "math/number_theory.hpp"

#include <iostream>
#include <vector>

int main() {
    auto [gcd, x, y] = m1une::math::extended_gcd(30, 18);
    std::cout << gcd << " " << x << " " << y << "\n"; // 6 -1 2

    auto solution = m1une::math::crt(
        std::vector<long long>{2, 3, 2},
        std::vector<long long>{3, 5, 7}
    );
    std::cout << solution.first << " " << solution.second << "\n";  // 23 105

    std::cout << m1une::math::inv_mod(3, 11) << "\n";  // 4
    std::cout << m1une::math::floor_sum(4, 10, 6, 3) << "\n";  // 3
}

Required by

Verified with

Code

#ifndef M1UNE_MATH_NUMBER_THEORY_HPP
#define M1UNE_MATH_NUMBER_THEORY_HPP 1

#include <cassert>
#include <cstdint>
#include <limits>
#include <tuple>
#include <utility>
#include <vector>

namespace m1une {
namespace math {

namespace internal {

inline long long safe_mod(long long x, long long mod) {
    x %= mod;
    if (x < 0) x += mod;
    return x;
}

inline unsigned __int128 floor_sum_unsigned(unsigned long long n, unsigned long long mod, unsigned long long a,
                                            unsigned long long b) {
    unsigned __int128 answer = 0;
    while (true) {
        if (a >= mod) {
            answer += static_cast<unsigned __int128>(n) * (n - 1) / 2 * (a / mod);
            a %= mod;
        }
        if (b >= mod) {
            answer += static_cast<unsigned __int128>(n) * (b / mod);
            b %= mod;
        }

        const unsigned __int128 y_max = static_cast<unsigned __int128>(a) * n + b;
        if (y_max < mod) break;
        n = static_cast<unsigned long long>(y_max / mod);
        b = static_cast<unsigned long long>(y_max % mod);
        unsigned long long tmp = mod;
        mod = a;
        a = tmp;
    }
    return answer;
}

}  // namespace internal

// Returns (g, x, y), where g = gcd(a, b) is nonnegative and
// a * x + b * y = g. Returns (0, 0, 0) when a = b = 0.
inline std::tuple<long long, long long, long long> extended_gcd(long long a,
                                                               long long b) {
    using i128 = __int128;
    if (a == 0 && b == 0) return {0, 0, 0};

    i128 old_remainder = a;
    i128 remainder = b;
    if (old_remainder < 0) old_remainder = -old_remainder;
    if (remainder < 0) remainder = -remainder;
    i128 old_x = 1;
    i128 x = 0;
    i128 old_y = 0;
    i128 y = 1;

    while (remainder != 0) {
        i128 quotient = old_remainder / remainder;

        i128 next = old_remainder - quotient * remainder;
        old_remainder = remainder;
        remainder = next;

        next = old_x - quotient * x;
        old_x = x;
        x = next;

        next = old_y - quotient * y;
        old_y = y;
        y = next;
    }

    if (a < 0) old_x = -old_x;
    if (b < 0) old_y = -old_y;

#ifndef NDEBUG
    const i128 minimum = std::numeric_limits<long long>::min();
    const i128 maximum = std::numeric_limits<long long>::max();
    assert(old_remainder <= maximum);
    assert(minimum <= old_x && old_x <= maximum);
    assert(minimum <= old_y && old_y <= maximum);
#endif
    return {static_cast<long long>(old_remainder), static_cast<long long>(old_x),
            static_cast<long long>(old_y)};
}

inline long long pow_mod(long long x, unsigned long long exponent, long long mod) {
    assert(mod >= 1);
    if (mod == 1) return 0;

    unsigned long long base = static_cast<unsigned long long>(internal::safe_mod(x, mod));
    unsigned long long result = 1;
    const unsigned long long unsigned_mod = static_cast<unsigned long long>(mod);
    while (exponent > 0) {
        if (exponent & 1) {
            result = static_cast<unsigned long long>(static_cast<unsigned __int128>(result) * base % unsigned_mod);
        }
        base = static_cast<unsigned long long>(static_cast<unsigned __int128>(base) * base % unsigned_mod);
        exponent >>= 1;
    }
    return static_cast<long long>(result);
}

// Returns gcd(a, mod) and x such that a * x is congruent to gcd(a, mod)
// modulo mod. The returned x is in [0, mod / gcd(a, mod)).
inline std::pair<long long, long long> inv_gcd(long long a, long long mod) {
    assert(mod >= 1);
    a = internal::safe_mod(a, mod);
    if (a == 0) return {mod, 0};

    long long s = mod;
    long long t = a;
    long long m0 = 0;
    long long m1 = 1;
    while (t > 0) {
        const long long quotient = s / t;
        s -= t * quotient;
        m0 -= m1 * quotient;

        long long tmp = s;
        s = t;
        t = tmp;
        tmp = m0;
        m0 = m1;
        m1 = tmp;
    }
    if (m0 < 0) m0 += mod / s;
    return {s, m0};
}

inline long long inv_mod(long long x, long long mod) {
    const auto result = inv_gcd(x, mod);
    assert(result.first == 1);
    return result.second;
}

// Returns the smallest nonnegative solution and the least common multiple of
// the moduli. Returns {0, 0} when the system is inconsistent.
inline std::pair<long long, long long> crt(const std::vector<long long>& remainders,
                                           const std::vector<long long>& moduli) {
    assert(remainders.size() == moduli.size());

    long long r0 = 0;
    long long m0 = 1;
    for (int i = 0; i < int(remainders.size()); i++) {
        assert(moduli[i] >= 1);
        long long r1 = internal::safe_mod(remainders[i], moduli[i]);
        long long m1 = moduli[i];

        if (m0 < m1) {
            long long tmp = r0;
            r0 = r1;
            r1 = tmp;
            tmp = m0;
            m0 = m1;
            m1 = tmp;
        }
        if (m0 % m1 == 0) {
            if (r0 % m1 != r1) return {0, 0};
            continue;
        }

        const auto inverse = inv_gcd(m0, m1);
        const long long gcd = inverse.first;
        const long long reduced_modulus = m1 / gcd;
        const __int128 difference = static_cast<__int128>(r1) - r0;
        if (difference % gcd != 0) return {0, 0};

        __int128 multiplier = difference / gcd % reduced_modulus;
        multiplier = multiplier * inverse.second % reduced_modulus;
        if (multiplier < 0) multiplier += reduced_modulus;

        const __int128 new_modulus = static_cast<__int128>(m0) * reduced_modulus;
        assert(new_modulus <= std::numeric_limits<long long>::max());
        __int128 new_remainder = static_cast<__int128>(r0) + multiplier * m0;
        new_remainder %= new_modulus;
        if (new_remainder < 0) new_remainder += new_modulus;
        r0 = static_cast<long long>(new_remainder);
        m0 = static_cast<long long>(new_modulus);
    }
    return {r0, m0};
}

// Returns sum_{i=0}^{n-1} floor((a * i + b) / mod).
inline long long floor_sum(long long n, long long mod, long long a, long long b) {
    assert(n >= 0);
    assert(mod >= 1);

    const long long normalized_a = internal::safe_mod(a, mod);
    const long long normalized_b = internal::safe_mod(b, mod);
    __int128 answer = (static_cast<__int128>(a) - normalized_a) / mod * n * (n - 1) / 2;
    answer += (static_cast<__int128>(b) - normalized_b) / mod * n;
    answer += internal::floor_sum_unsigned(static_cast<unsigned long long>(n), static_cast<unsigned long long>(mod),
                                           static_cast<unsigned long long>(normalized_a),
                                           static_cast<unsigned long long>(normalized_b));

    assert(answer >= std::numeric_limits<long long>::min());
    assert(answer <= std::numeric_limits<long long>::max());
    return static_cast<long long>(answer);
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_NUMBER_THEORY_HPP

#line 1 "math/number_theory.hpp"



#include <cassert>
#include <cstdint>
#include <limits>
#include <tuple>
#include <utility>
#include <vector>

namespace m1une {
namespace math {

namespace internal {

inline long long safe_mod(long long x, long long mod) {
    x %= mod;
    if (x < 0) x += mod;
    return x;
}

inline unsigned __int128 floor_sum_unsigned(unsigned long long n, unsigned long long mod, unsigned long long a,
                                            unsigned long long b) {
    unsigned __int128 answer = 0;
    while (true) {
        if (a >= mod) {
            answer += static_cast<unsigned __int128>(n) * (n - 1) / 2 * (a / mod);
            a %= mod;
        }
        if (b >= mod) {
            answer += static_cast<unsigned __int128>(n) * (b / mod);
            b %= mod;
        }

        const unsigned __int128 y_max = static_cast<unsigned __int128>(a) * n + b;
        if (y_max < mod) break;
        n = static_cast<unsigned long long>(y_max / mod);
        b = static_cast<unsigned long long>(y_max % mod);
        unsigned long long tmp = mod;
        mod = a;
        a = tmp;
    }
    return answer;
}

}  // namespace internal

// Returns (g, x, y), where g = gcd(a, b) is nonnegative and
// a * x + b * y = g. Returns (0, 0, 0) when a = b = 0.
inline std::tuple<long long, long long, long long> extended_gcd(long long a,
                                                               long long b) {
    using i128 = __int128;
    if (a == 0 && b == 0) return {0, 0, 0};

    i128 old_remainder = a;
    i128 remainder = b;
    if (old_remainder < 0) old_remainder = -old_remainder;
    if (remainder < 0) remainder = -remainder;
    i128 old_x = 1;
    i128 x = 0;
    i128 old_y = 0;
    i128 y = 1;

    while (remainder != 0) {
        i128 quotient = old_remainder / remainder;

        i128 next = old_remainder - quotient * remainder;
        old_remainder = remainder;
        remainder = next;

        next = old_x - quotient * x;
        old_x = x;
        x = next;

        next = old_y - quotient * y;
        old_y = y;
        y = next;
    }

    if (a < 0) old_x = -old_x;
    if (b < 0) old_y = -old_y;

#ifndef NDEBUG
    const i128 minimum = std::numeric_limits<long long>::min();
    const i128 maximum = std::numeric_limits<long long>::max();
    assert(old_remainder <= maximum);
    assert(minimum <= old_x && old_x <= maximum);
    assert(minimum <= old_y && old_y <= maximum);
#endif
    return {static_cast<long long>(old_remainder), static_cast<long long>(old_x),
            static_cast<long long>(old_y)};
}

inline long long pow_mod(long long x, unsigned long long exponent, long long mod) {
    assert(mod >= 1);
    if (mod == 1) return 0;

    unsigned long long base = static_cast<unsigned long long>(internal::safe_mod(x, mod));
    unsigned long long result = 1;
    const unsigned long long unsigned_mod = static_cast<unsigned long long>(mod);
    while (exponent > 0) {
        if (exponent & 1) {
            result = static_cast<unsigned long long>(static_cast<unsigned __int128>(result) * base % unsigned_mod);
        }
        base = static_cast<unsigned long long>(static_cast<unsigned __int128>(base) * base % unsigned_mod);
        exponent >>= 1;
    }
    return static_cast<long long>(result);
}

// Returns gcd(a, mod) and x such that a * x is congruent to gcd(a, mod)
// modulo mod. The returned x is in [0, mod / gcd(a, mod)).
inline std::pair<long long, long long> inv_gcd(long long a, long long mod) {
    assert(mod >= 1);
    a = internal::safe_mod(a, mod);
    if (a == 0) return {mod, 0};

    long long s = mod;
    long long t = a;
    long long m0 = 0;
    long long m1 = 1;
    while (t > 0) {
        const long long quotient = s / t;
        s -= t * quotient;
        m0 -= m1 * quotient;

        long long tmp = s;
        s = t;
        t = tmp;
        tmp = m0;
        m0 = m1;
        m1 = tmp;
    }
    if (m0 < 0) m0 += mod / s;
    return {s, m0};
}

inline long long inv_mod(long long x, long long mod) {
    const auto result = inv_gcd(x, mod);
    assert(result.first == 1);
    return result.second;
}

// Returns the smallest nonnegative solution and the least common multiple of
// the moduli. Returns {0, 0} when the system is inconsistent.
inline std::pair<long long, long long> crt(const std::vector<long long>& remainders,
                                           const std::vector<long long>& moduli) {
    assert(remainders.size() == moduli.size());

    long long r0 = 0;
    long long m0 = 1;
    for (int i = 0; i < int(remainders.size()); i++) {
        assert(moduli[i] >= 1);
        long long r1 = internal::safe_mod(remainders[i], moduli[i]);
        long long m1 = moduli[i];

        if (m0 < m1) {
            long long tmp = r0;
            r0 = r1;
            r1 = tmp;
            tmp = m0;
            m0 = m1;
            m1 = tmp;
        }
        if (m0 % m1 == 0) {
            if (r0 % m1 != r1) return {0, 0};
            continue;
        }

        const auto inverse = inv_gcd(m0, m1);
        const long long gcd = inverse.first;
        const long long reduced_modulus = m1 / gcd;
        const __int128 difference = static_cast<__int128>(r1) - r0;
        if (difference % gcd != 0) return {0, 0};

        __int128 multiplier = difference / gcd % reduced_modulus;
        multiplier = multiplier * inverse.second % reduced_modulus;
        if (multiplier < 0) multiplier += reduced_modulus;

        const __int128 new_modulus = static_cast<__int128>(m0) * reduced_modulus;
        assert(new_modulus <= std::numeric_limits<long long>::max());
        __int128 new_remainder = static_cast<__int128>(r0) + multiplier * m0;
        new_remainder %= new_modulus;
        if (new_remainder < 0) new_remainder += new_modulus;
        r0 = static_cast<long long>(new_remainder);
        m0 = static_cast<long long>(new_modulus);
    }
    return {r0, m0};
}

// Returns sum_{i=0}^{n-1} floor((a * i + b) / mod).
inline long long floor_sum(long long n, long long mod, long long a, long long b) {
    assert(n >= 0);
    assert(mod >= 1);

    const long long normalized_a = internal::safe_mod(a, mod);
    const long long normalized_b = internal::safe_mod(b, mod);
    __int128 answer = (static_cast<__int128>(a) - normalized_a) / mod * n * (n - 1) / 2;
    answer += (static_cast<__int128>(b) - normalized_b) / mod * n;
    answer += internal::floor_sum_unsigned(static_cast<unsigned long long>(n), static_cast<unsigned long long>(mod),
                                           static_cast<unsigned long long>(normalized_a),
                                           static_cast<unsigned long long>(normalized_b));

    assert(answer >= std::numeric_limits<long long>::min());
    assert(answer <= std::numeric_limits<long long>::max());
    return static_cast<long long>(answer);
}

}  // namespace math
}  // namespace m1une

m1une's library

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

Number Theory
(math/number_theory.hpp)

Overview

Extended Greatest Common Divisor

Congruences and Modular Arithmetic

Modular Power

Modular Inverse

Chinese Remainder Theorem

Floor Sum

API

Example

Required by

Verified with

Code

Number Theory (math/number_theory.hpp)

Overview

Extended Greatest Common Divisor

Congruences and Modular Arithmetic

Modular Power

Modular Inverse

Chinese Remainder Theorem

Floor Sum

API

Example

Required by

Verified with

Code

Number Theory
(math/number_theory.hpp)