m1une's library

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

View on GitHub

:heavy_check_mark: Gaussian Integer
(math/gaussian_integer.hpp)

Overview

GaussianInteger<T> represents an exact number a + b i with integral real and imaginary parts. Gaussian integers form a Euclidean domain under the norm

\[N(a + b i) = a^2 + b^2.\]

The type provides ordinary arithmetic, conjugation, exact nearest-quotient Euclidean division, canonical associates, divisibility tests, GCD, and extended GCD. All calculations are integral; division never uses floating point.

#include "math/gaussian_integer.hpp"

Template Requirements

T must be a standard signed integral type no wider than long long. The default is long long. Coordinates returned by arithmetic and extended GCD must remain representable by T; debug builds assert this condition.

Products, norms, and division numerators use exact 128-bit intermediates. The public norm_type is __uint128_t, so the norm of every pair of 64-bit coordinates is representable.

Representation and Fields

Member Meaning
value_type Alias for T.
norm_type Unsigned 128-bit norm type.
real Real coordinate a.
imag Imaginary coordinate b.

Equality compares both coordinates exactly.

Euclidean Division

For a nonzero divisor d, z.divmod(d) returns (q, r) satisfying

\[z = qd + r, \qquad N(r) < N(d).\]

The real and imaginary coordinates of the exact complex quotient are rounded to their nearest integers. Exact half-integer ties are rounded away from zero; any consistent nearest choice gives a valid Euclidean remainder.

Consequently, / is Euclidean nearest-quotient division, not division in the complex field. % returns its matching remainder. Division by zero is invalid and asserted.

Canonical Associates

The units are 1, -1, i, and -i. Multiplying by a unit does not change divisibility, so a nonzero Gaussian GCD has four associates.

normalized() chooses the unique associate with positive real part and nonnegative imaginary part. Zero remains zero. gaussian_gcd and extended_gaussian_gcd return this canonical associate, making their output deterministic.

Interface

All members and free functions are in m1une::math.

Method Exact signature Description Complexity
Default constructor GaussianInteger() Constructs zero. $O(1)$
Real constructor GaussianInteger(T real) Constructs a real Gaussian integer. $O(1)$
Coordinate constructor GaussianInteger(T real, T imag) Constructs real + imag * i. $O(1)$
is_zero bool is_zero() const Tests equality with zero. $O(1)$
is_unit bool is_unit() const Tests whether the norm is one. $O(1)$
norm norm_type norm() const Returns real * real + imag * imag. $O(1)$
conjugate GaussianInteger conjugate() const Returns real - imag * i. $O(1)$
normalizing_unit GaussianInteger normalizing_unit() const Returns the unit that moves the value to its canonical associate. $O(1)$
normalized GaussianInteger normalized() const Returns the canonical associate. $O(1)$
divmod std::pair<GaussianInteger, GaussianInteger> divmod(const GaussianInteger& divisor) const Returns the Euclidean quotient and remainder. $O(1)$
Unary plus GaussianInteger operator+() const Returns the unchanged value. $O(1)$
Unary minus GaussianInteger operator-() const Negates both coordinates. $O(1)$
Compound addition GaussianInteger& operator+=(const GaussianInteger& other) Adds other. $O(1)$
Compound subtraction GaussianInteger& operator-=(const GaussianInteger& other) Subtracts other. $O(1)$
Compound multiplication GaussianInteger& operator*=(const GaussianInteger& other) Multiplies by other. $O(1)$
Compound division GaussianInteger& operator/=(const GaussianInteger& other) Replaces the value by its Euclidean quotient by other. $O(1)$
Compound remainder GaussianInteger& operator%=(const GaussianInteger& other) Replaces the value by its Euclidean remainder modulo other. $O(1)$
Binary addition friend GaussianInteger operator+(GaussianInteger left, const GaussianInteger& right) Returns left + right. $O(1)$
Binary subtraction friend GaussianInteger operator-(GaussianInteger left, const GaussianInteger& right) Returns left - right. $O(1)$
Binary multiplication friend GaussianInteger operator*(GaussianInteger left, const GaussianInteger& right) Returns left * right. $O(1)$
Binary division friend GaussianInteger operator/(GaussianInteger left, const GaussianInteger& right) Returns the Euclidean quotient. $O(1)$
Binary remainder friend GaussianInteger operator%(GaussianInteger left, const GaussianInteger& right) Returns the Euclidean remainder. $O(1)$
Equality friend bool operator==(const GaussianInteger& first, const GaussianInteger& second) Compares both coordinates. $O(1)$
Free function Exact signature Description Complexity
gaussian_divides bool gaussian_divides(const GaussianInteger<T>& divisor, const GaussianInteger<T>& value) Tests divisor dividing value; zero divides only zero. $O(1)$
gaussian_associates bool gaussian_associates(const GaussianInteger<T>& first, const GaussianInteger<T>& second) Tests equality up to multiplication by a unit. $O(1)$
gaussian_gcd GaussianInteger<T> gaussian_gcd(GaussianInteger<T> first, GaussianInteger<T> second) Returns the canonical greatest common divisor. $O(\log(\min(N(first), N(second)) + 1))$
extended_gaussian_gcd std::tuple<GaussianInteger<T>, GaussianInteger<T>, GaussianInteger<T>> extended_gaussian_gcd(GaussianInteger<T> first, GaussianInteger<T> second) Returns (g, x, y) with first * x + second * y == g, where g is canonical. $O(\log(\min(N(first), N(second)) + 1))$

The zero pair has GCD zero. For gaussian_divides, the convention that zero divides zero follows directly from the existence of a Gaussian integer z with 0 = 0 * z.

Example

#include "math/gaussian_integer.hpp"

#include <iostream>

int main() {
    using G = m1une::math::GaussianInteger<long long>;

    G first(5, 7);
    G second(2, -1);
    auto [quotient, remainder] = first.divmod(second);

    G gcd = m1une::math::gaussian_gcd(first, second);
    std::cout << quotient.real << ' ' << quotient.imag << '\n';
    std::cout << remainder.real << ' ' << remainder.imag << '\n';
    std::cout << gcd.real << ' ' << gcd.imag << '\n';
}

Required by

Verified with

Code

#ifndef M1UNE_MATH_GAUSSIAN_INTEGER_HPP
#define M1UNE_MATH_GAUSSIAN_INTEGER_HPP 1

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

namespace m1une {
namespace math {

template <std::signed_integral T = long long>
struct GaussianInteger {
    static_assert(sizeof(T) <= sizeof(long long));

    using value_type = T;
    using norm_type = __uint128_t;

    T real;
    T imag;

   private:
    using wide_type = __int128_t;

    struct SignedMagnitude {
        norm_type magnitude;
        bool negative;
    };

    static constexpr norm_type magnitude(T value) {
        const wide_type wide = value;
        if (wide < 0) {
            return static_cast<norm_type>(-(wide + 1)) + 1;
        }
        return static_cast<norm_type>(wide);
    }

    static constexpr SignedMagnitude signed_value(T value) {
        return SignedMagnitude{magnitude(value), value < 0};
    }

    static constexpr SignedMagnitude negate(SignedMagnitude value) {
        if (value.magnitude != 0) value.negative = !value.negative;
        return value;
    }

    static constexpr SignedMagnitude add_signed(
        SignedMagnitude first,
        SignedMagnitude second
    ) {
        if (first.negative == second.negative) {
            return SignedMagnitude{
                first.magnitude + second.magnitude,
                first.negative
            };
        }
        if (first.magnitude < second.magnitude) {
            return SignedMagnitude{
                second.magnitude - first.magnitude,
                second.negative
            };
        }
        return SignedMagnitude{
            first.magnitude - second.magnitude,
            first.magnitude == second.magnitude ? false : first.negative
        };
    }

    static constexpr SignedMagnitude product(T first, T second) {
        const norm_type result = magnitude(first) * magnitude(second);
        return SignedMagnitude{
            result,
            result != 0 && ((first < 0) != (second < 0))
        };
    }

    static constexpr T narrow(SignedMagnitude value) {
        const norm_type maximum =
            static_cast<norm_type>(std::numeric_limits<T>::max());
        if (!value.negative) {
            assert(value.magnitude <= maximum);
            return static_cast<T>(value.magnitude);
        }

        assert(value.magnitude <= maximum + 1);
        if (value.magnitude == maximum + 1) {
            return std::numeric_limits<T>::min();
        }
        return static_cast<T>(-static_cast<wide_type>(value.magnitude));
    }

    static constexpr std::pair<SignedMagnitude, SignedMagnitude>
    product_components(
        const GaussianInteger& first,
        const GaussianInteger& second
    ) {
        const SignedMagnitude product_real = add_signed(
            product(first.real, second.real),
            negate(product(first.imag, second.imag))
        );
        const SignedMagnitude product_imag = add_signed(
            product(first.real, second.imag),
            product(first.imag, second.real)
        );
        return {product_real, product_imag};
    }

    static constexpr T round_ratio(
        SignedMagnitude numerator,
        norm_type denominator
    ) {
        assert(denominator != 0);
        norm_type quotient = numerator.magnitude / denominator;
        const norm_type remainder = numerator.magnitude % denominator;
        if (remainder + remainder >= denominator) quotient++;
        return narrow(SignedMagnitude{quotient, numerator.negative});
    }

   public:
    constexpr GaussianInteger() : real(0), imag(0) {}

    constexpr GaussianInteger(T real_) : real(real_), imag(0) {}

    constexpr GaussianInteger(T real_, T imag_)
        : real(real_), imag(imag_) {}

    constexpr bool is_zero() const {
        return real == 0 && imag == 0;
    }

    constexpr bool is_unit() const {
        return norm() == 1;
    }

    constexpr norm_type norm() const {
        const norm_type real_magnitude = magnitude(real);
        const norm_type imag_magnitude = magnitude(imag);
        return real_magnitude * real_magnitude +
               imag_magnitude * imag_magnitude;
    }

    constexpr GaussianInteger conjugate() const {
        return GaussianInteger(real, narrow(negate(signed_value(imag))));
    }

    constexpr GaussianInteger normalizing_unit() const {
        if (is_zero() || (0 < real && 0 <= imag)) {
            return GaussianInteger(1, 0);
        }
        if (0 < real) return GaussianInteger(0, 1);
        if (real < 0 && imag <= 0) return GaussianInteger(-1, 0);
        if (real < 0) return GaussianInteger(0, -1);
        if (0 < imag) return GaussianInteger(0, -1);
        return GaussianInteger(0, 1);
    }

    constexpr GaussianInteger normalized() const {
        return normalizing_unit() * *this;
    }

    constexpr std::pair<GaussianInteger, GaussianInteger> divmod(
        const GaussianInteger& divisor
    ) const {
        const norm_type divisor_norm = divisor.norm();
        assert(divisor_norm != 0);

        const SignedMagnitude numerator_real = add_signed(
            product(real, divisor.real),
            product(imag, divisor.imag)
        );
        const SignedMagnitude numerator_imag = add_signed(
            product(imag, divisor.real),
            negate(product(real, divisor.imag))
        );
        const GaussianInteger quotient(
            round_ratio(numerator_real, divisor_norm),
            round_ratio(numerator_imag, divisor_norm)
        );

        const auto product = product_components(quotient, divisor);
        const GaussianInteger remainder(
            narrow(add_signed(signed_value(real), negate(product.first))),
            narrow(add_signed(signed_value(imag), negate(product.second)))
        );
        assert(remainder.norm() < divisor_norm);
        return {quotient, remainder};
    }

    constexpr GaussianInteger operator+() const {
        return *this;
    }

    constexpr GaussianInteger operator-() const {
        return GaussianInteger(
            narrow(negate(signed_value(real))),
            narrow(negate(signed_value(imag)))
        );
    }

    constexpr GaussianInteger& operator+=(const GaussianInteger& other) {
        real = narrow(add_signed(signed_value(real), signed_value(other.real)));
        imag = narrow(add_signed(signed_value(imag), signed_value(other.imag)));
        return *this;
    }

    constexpr GaussianInteger& operator-=(const GaussianInteger& other) {
        real = narrow(add_signed(
            signed_value(real),
            negate(signed_value(other.real))
        ));
        imag = narrow(add_signed(
            signed_value(imag),
            negate(signed_value(other.imag))
        ));
        return *this;
    }

    constexpr GaussianInteger& operator*=(const GaussianInteger& other) {
        const auto result = product_components(*this, other);
        real = narrow(result.first);
        imag = narrow(result.second);
        return *this;
    }

    constexpr GaussianInteger& operator/=(const GaussianInteger& other) {
        *this = divmod(other).first;
        return *this;
    }

    constexpr GaussianInteger& operator%=(const GaussianInteger& other) {
        *this = divmod(other).second;
        return *this;
    }

    friend constexpr GaussianInteger operator+(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left += right;
    }

    friend constexpr GaussianInteger operator-(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left -= right;
    }

    friend constexpr GaussianInteger operator*(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left *= right;
    }

    friend constexpr GaussianInteger operator/(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left /= right;
    }

    friend constexpr GaussianInteger operator%(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left %= right;
    }

    friend constexpr bool operator==(
        const GaussianInteger& first,
        const GaussianInteger& second
    ) = default;
};

template <std::signed_integral T>
constexpr bool gaussian_divides(
    const GaussianInteger<T>& divisor,
    const GaussianInteger<T>& value
) {
    if (divisor.is_zero()) return value.is_zero();
    return (value % divisor).is_zero();
}

template <std::signed_integral T>
constexpr bool gaussian_associates(
    const GaussianInteger<T>& first,
    const GaussianInteger<T>& second
) {
    return first.normalized() == second.normalized();
}

template <std::signed_integral T>
constexpr GaussianInteger<T> gaussian_gcd(
    GaussianInteger<T> first,
    GaussianInteger<T> second
) {
    while (!second.is_zero()) {
        first %= second;
        std::swap(first, second);
    }
    return first.normalized();
}

template <std::signed_integral T>
constexpr std::tuple<
    GaussianInteger<T>,
    GaussianInteger<T>,
    GaussianInteger<T>
> extended_gaussian_gcd(
    GaussianInteger<T> first,
    GaussianInteger<T> second
) {
    using G = GaussianInteger<T>;
    G old_remainder = first;
    G remainder = second;
    G old_first_coefficient(1);
    G first_coefficient(0);
    G old_second_coefficient(0);
    G second_coefficient(1);

    while (!remainder.is_zero()) {
        const G quotient = old_remainder / remainder;

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

        next = old_first_coefficient - quotient * first_coefficient;
        old_first_coefficient = first_coefficient;
        first_coefficient = next;

        next = old_second_coefficient - quotient * second_coefficient;
        old_second_coefficient = second_coefficient;
        second_coefficient = next;
    }

    const G unit = old_remainder.normalizing_unit();
    return {
        unit * old_remainder,
        unit * old_first_coefficient,
        unit * old_second_coefficient
    };
}

}  // namespace math
}  // namespace m1une

#endif  // M1UNE_MATH_GAUSSIAN_INTEGER_HPP
#line 1 "math/gaussian_integer.hpp"



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

namespace m1une {
namespace math {

template <std::signed_integral T = long long>
struct GaussianInteger {
    static_assert(sizeof(T) <= sizeof(long long));

    using value_type = T;
    using norm_type = __uint128_t;

    T real;
    T imag;

   private:
    using wide_type = __int128_t;

    struct SignedMagnitude {
        norm_type magnitude;
        bool negative;
    };

    static constexpr norm_type magnitude(T value) {
        const wide_type wide = value;
        if (wide < 0) {
            return static_cast<norm_type>(-(wide + 1)) + 1;
        }
        return static_cast<norm_type>(wide);
    }

    static constexpr SignedMagnitude signed_value(T value) {
        return SignedMagnitude{magnitude(value), value < 0};
    }

    static constexpr SignedMagnitude negate(SignedMagnitude value) {
        if (value.magnitude != 0) value.negative = !value.negative;
        return value;
    }

    static constexpr SignedMagnitude add_signed(
        SignedMagnitude first,
        SignedMagnitude second
    ) {
        if (first.negative == second.negative) {
            return SignedMagnitude{
                first.magnitude + second.magnitude,
                first.negative
            };
        }
        if (first.magnitude < second.magnitude) {
            return SignedMagnitude{
                second.magnitude - first.magnitude,
                second.negative
            };
        }
        return SignedMagnitude{
            first.magnitude - second.magnitude,
            first.magnitude == second.magnitude ? false : first.negative
        };
    }

    static constexpr SignedMagnitude product(T first, T second) {
        const norm_type result = magnitude(first) * magnitude(second);
        return SignedMagnitude{
            result,
            result != 0 && ((first < 0) != (second < 0))
        };
    }

    static constexpr T narrow(SignedMagnitude value) {
        const norm_type maximum =
            static_cast<norm_type>(std::numeric_limits<T>::max());
        if (!value.negative) {
            assert(value.magnitude <= maximum);
            return static_cast<T>(value.magnitude);
        }

        assert(value.magnitude <= maximum + 1);
        if (value.magnitude == maximum + 1) {
            return std::numeric_limits<T>::min();
        }
        return static_cast<T>(-static_cast<wide_type>(value.magnitude));
    }

    static constexpr std::pair<SignedMagnitude, SignedMagnitude>
    product_components(
        const GaussianInteger& first,
        const GaussianInteger& second
    ) {
        const SignedMagnitude product_real = add_signed(
            product(first.real, second.real),
            negate(product(first.imag, second.imag))
        );
        const SignedMagnitude product_imag = add_signed(
            product(first.real, second.imag),
            product(first.imag, second.real)
        );
        return {product_real, product_imag};
    }

    static constexpr T round_ratio(
        SignedMagnitude numerator,
        norm_type denominator
    ) {
        assert(denominator != 0);
        norm_type quotient = numerator.magnitude / denominator;
        const norm_type remainder = numerator.magnitude % denominator;
        if (remainder + remainder >= denominator) quotient++;
        return narrow(SignedMagnitude{quotient, numerator.negative});
    }

   public:
    constexpr GaussianInteger() : real(0), imag(0) {}

    constexpr GaussianInteger(T real_) : real(real_), imag(0) {}

    constexpr GaussianInteger(T real_, T imag_)
        : real(real_), imag(imag_) {}

    constexpr bool is_zero() const {
        return real == 0 && imag == 0;
    }

    constexpr bool is_unit() const {
        return norm() == 1;
    }

    constexpr norm_type norm() const {
        const norm_type real_magnitude = magnitude(real);
        const norm_type imag_magnitude = magnitude(imag);
        return real_magnitude * real_magnitude +
               imag_magnitude * imag_magnitude;
    }

    constexpr GaussianInteger conjugate() const {
        return GaussianInteger(real, narrow(negate(signed_value(imag))));
    }

    constexpr GaussianInteger normalizing_unit() const {
        if (is_zero() || (0 < real && 0 <= imag)) {
            return GaussianInteger(1, 0);
        }
        if (0 < real) return GaussianInteger(0, 1);
        if (real < 0 && imag <= 0) return GaussianInteger(-1, 0);
        if (real < 0) return GaussianInteger(0, -1);
        if (0 < imag) return GaussianInteger(0, -1);
        return GaussianInteger(0, 1);
    }

    constexpr GaussianInteger normalized() const {
        return normalizing_unit() * *this;
    }

    constexpr std::pair<GaussianInteger, GaussianInteger> divmod(
        const GaussianInteger& divisor
    ) const {
        const norm_type divisor_norm = divisor.norm();
        assert(divisor_norm != 0);

        const SignedMagnitude numerator_real = add_signed(
            product(real, divisor.real),
            product(imag, divisor.imag)
        );
        const SignedMagnitude numerator_imag = add_signed(
            product(imag, divisor.real),
            negate(product(real, divisor.imag))
        );
        const GaussianInteger quotient(
            round_ratio(numerator_real, divisor_norm),
            round_ratio(numerator_imag, divisor_norm)
        );

        const auto product = product_components(quotient, divisor);
        const GaussianInteger remainder(
            narrow(add_signed(signed_value(real), negate(product.first))),
            narrow(add_signed(signed_value(imag), negate(product.second)))
        );
        assert(remainder.norm() < divisor_norm);
        return {quotient, remainder};
    }

    constexpr GaussianInteger operator+() const {
        return *this;
    }

    constexpr GaussianInteger operator-() const {
        return GaussianInteger(
            narrow(negate(signed_value(real))),
            narrow(negate(signed_value(imag)))
        );
    }

    constexpr GaussianInteger& operator+=(const GaussianInteger& other) {
        real = narrow(add_signed(signed_value(real), signed_value(other.real)));
        imag = narrow(add_signed(signed_value(imag), signed_value(other.imag)));
        return *this;
    }

    constexpr GaussianInteger& operator-=(const GaussianInteger& other) {
        real = narrow(add_signed(
            signed_value(real),
            negate(signed_value(other.real))
        ));
        imag = narrow(add_signed(
            signed_value(imag),
            negate(signed_value(other.imag))
        ));
        return *this;
    }

    constexpr GaussianInteger& operator*=(const GaussianInteger& other) {
        const auto result = product_components(*this, other);
        real = narrow(result.first);
        imag = narrow(result.second);
        return *this;
    }

    constexpr GaussianInteger& operator/=(const GaussianInteger& other) {
        *this = divmod(other).first;
        return *this;
    }

    constexpr GaussianInteger& operator%=(const GaussianInteger& other) {
        *this = divmod(other).second;
        return *this;
    }

    friend constexpr GaussianInteger operator+(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left += right;
    }

    friend constexpr GaussianInteger operator-(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left -= right;
    }

    friend constexpr GaussianInteger operator*(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left *= right;
    }

    friend constexpr GaussianInteger operator/(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left /= right;
    }

    friend constexpr GaussianInteger operator%(
        GaussianInteger left,
        const GaussianInteger& right
    ) {
        return left %= right;
    }

    friend constexpr bool operator==(
        const GaussianInteger& first,
        const GaussianInteger& second
    ) = default;
};

template <std::signed_integral T>
constexpr bool gaussian_divides(
    const GaussianInteger<T>& divisor,
    const GaussianInteger<T>& value
) {
    if (divisor.is_zero()) return value.is_zero();
    return (value % divisor).is_zero();
}

template <std::signed_integral T>
constexpr bool gaussian_associates(
    const GaussianInteger<T>& first,
    const GaussianInteger<T>& second
) {
    return first.normalized() == second.normalized();
}

template <std::signed_integral T>
constexpr GaussianInteger<T> gaussian_gcd(
    GaussianInteger<T> first,
    GaussianInteger<T> second
) {
    while (!second.is_zero()) {
        first %= second;
        std::swap(first, second);
    }
    return first.normalized();
}

template <std::signed_integral T>
constexpr std::tuple<
    GaussianInteger<T>,
    GaussianInteger<T>,
    GaussianInteger<T>
> extended_gaussian_gcd(
    GaussianInteger<T> first,
    GaussianInteger<T> second
) {
    using G = GaussianInteger<T>;
    G old_remainder = first;
    G remainder = second;
    G old_first_coefficient(1);
    G first_coefficient(0);
    G old_second_coefficient(0);
    G second_coefficient(1);

    while (!remainder.is_zero()) {
        const G quotient = old_remainder / remainder;

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

        next = old_first_coefficient - quotient * first_coefficient;
        old_first_coefficient = first_coefficient;
        first_coefficient = next;

        next = old_second_coefficient - quotient * second_coefficient;
        old_second_coefficient = second_coefficient;
        second_coefficient = next;
    }

    const G unit = old_remainder.normalizing_unit();
    return {
        unit * old_remainder,
        unit * old_first_coefficient,
        unit * old_second_coefficient
    };
}

}  // namespace math
}  // namespace m1une
Back to top page