gaussian_integer.hpp¶
Euclidean division in Z[i], returning quotient and remainder with norm(remainder) < norm(divisor).
Verified by gcd_of_gaussian_integers.
\[
\displaystyle \\frac{a+bi}{c+di}
\]
Implementation¶
#ifndef NOYA_GAUSSIAN_INTEGER_HPP
#define NOYA_GAUSSIAN_INTEGER_HPP 1
/// @complexity Time: O(log max(norm(a),norm(b))) Euclidean steps.
/// Space: O(1).
#include <cassert>
#include <cstdint>
#include <cstdlib>
#include <utility>
namespace noya {
struct gaussian_integer {
std::int64_t real = 0;
std::int64_t imaginary = 0;
friend bool operator==(const gaussian_integer &,
const gaussian_integer &) = default;
gaussian_integer &operator+=(const gaussian_integer &other) {
real += other.real;
imaginary += other.imaginary;
return *this;
}
gaussian_integer &operator-=(const gaussian_integer &other) {
real -= other.real;
imaginary -= other.imaginary;
return *this;
}
gaussian_integer &operator*=(const gaussian_integer &other) {
__int128 next_real =
__int128(real) * other.real - __int128(imaginary) * other.imaginary;
__int128 next_imaginary =
__int128(real) * other.imaginary + __int128(imaginary) * other.real;
real = std::int64_t(next_real);
imaginary = std::int64_t(next_imaginary);
return *this;
}
friend gaussian_integer operator+(gaussian_integer first,
const gaussian_integer &second) {
return first += second;
}
friend gaussian_integer operator-(gaussian_integer first,
const gaussian_integer &second) {
return first -= second;
}
friend gaussian_integer operator*(gaussian_integer first,
const gaussian_integer &second) {
return first *= second;
}
friend gaussian_integer operator-(gaussian_integer value) {
return {-value.real, -value.imaginary};
}
__int128 norm() const {
return __int128(real) * real + __int128(imaginary) * imaginary;
}
};
inline std::int64_t gaussian_round_div(__int128 numerator,
__int128 denominator) {
assert(denominator > 0);
__int128 quotient = numerator / denominator;
__int128 remainder = numerator % denominator;
__int128 absolute_remainder = remainder < 0 ? -remainder : remainder;
if (2 * absolute_remainder >= denominator) {
quotient += numerator >= 0 ? 1 : -1;
}
return std::int64_t(quotient);
}
/// @brief Euclidean division in Z[i], returning quotient and remainder with
/// norm(remainder) < norm(divisor).
inline std::pair<gaussian_integer, gaussian_integer>
gaussian_divmod(gaussian_integer dividend, gaussian_integer divisor) {
__int128 denominator = divisor.norm();
assert(denominator != 0);
__int128 real_numerator = __int128(dividend.real) * divisor.real +
__int128(dividend.imaginary) * divisor.imaginary;
__int128 imaginary_numerator = __int128(dividend.imaginary) * divisor.real -
__int128(dividend.real) * divisor.imaginary;
gaussian_integer quotient{
gaussian_round_div(real_numerator, denominator),
gaussian_round_div(imaginary_numerator, denominator)};
gaussian_integer remainder = dividend - quotient * divisor;
assert(remainder.norm() < divisor.norm());
return {quotient, remainder};
}
/// @brief Gaussian-integer gcd, normalized to positive real part, or positive
/// imaginary part when real is zero; associates differ only by a unit.
inline gaussian_integer gaussian_gcd(gaussian_integer first,
gaussian_integer second) {
while (second != gaussian_integer{}) {
auto [quotient, remainder] = gaussian_divmod(first, second);
(void)quotient;
first = second;
second = remainder;
}
if (first.real < 0 || (first.real == 0 && first.imaginary < 0)) {
first = -first;
}
return first;
}
} // namespace noya
#endif // NOYA_GAUSSIAN_INTEGER_HPP
#include <cassert>
#include <cstdint>
#include <cstdlib>
#include <utility>
/// @complexity Time: O(log max(norm(a),norm(b))) Euclidean steps.
/// Space: O(1).
namespace noya {
struct gaussian_integer {
std::int64_t real = 0;
std::int64_t imaginary = 0;
friend bool operator==(const gaussian_integer &,
const gaussian_integer &) = default;
gaussian_integer &operator+=(const gaussian_integer &other) {
real += other.real;
imaginary += other.imaginary;
return *this;
}
gaussian_integer &operator-=(const gaussian_integer &other) {
real -= other.real;
imaginary -= other.imaginary;
return *this;
}
gaussian_integer &operator*=(const gaussian_integer &other) {
__int128 next_real =
__int128(real) * other.real - __int128(imaginary) * other.imaginary;
__int128 next_imaginary =
__int128(real) * other.imaginary + __int128(imaginary) * other.real;
real = std::int64_t(next_real);
imaginary = std::int64_t(next_imaginary);
return *this;
}
friend gaussian_integer operator+(gaussian_integer first,
const gaussian_integer &second) {
return first += second;
}
friend gaussian_integer operator-(gaussian_integer first,
const gaussian_integer &second) {
return first -= second;
}
friend gaussian_integer operator*(gaussian_integer first,
const gaussian_integer &second) {
return first *= second;
}
friend gaussian_integer operator-(gaussian_integer value) {
return {-value.real, -value.imaginary};
}
__int128 norm() const {
return __int128(real) * real + __int128(imaginary) * imaginary;
}
};
inline std::int64_t gaussian_round_div(__int128 numerator,
__int128 denominator) {
assert(denominator > 0);
__int128 quotient = numerator / denominator;
__int128 remainder = numerator % denominator;
__int128 absolute_remainder = remainder < 0 ? -remainder : remainder;
if (2 * absolute_remainder >= denominator) {
quotient += numerator >= 0 ? 1 : -1;
}
return std::int64_t(quotient);
}
/// @brief Euclidean division in Z[i], returning quotient and remainder with
/// norm(remainder) < norm(divisor).
inline std::pair<gaussian_integer, gaussian_integer>
gaussian_divmod(gaussian_integer dividend, gaussian_integer divisor) {
__int128 denominator = divisor.norm();
assert(denominator != 0);
__int128 real_numerator = __int128(dividend.real) * divisor.real +
__int128(dividend.imaginary) * divisor.imaginary;
__int128 imaginary_numerator = __int128(dividend.imaginary) * divisor.real -
__int128(dividend.real) * divisor.imaginary;
gaussian_integer quotient{
gaussian_round_div(real_numerator, denominator),
gaussian_round_div(imaginary_numerator, denominator)};
gaussian_integer remainder = dividend - quotient * divisor;
assert(remainder.norm() < divisor.norm());
return {quotient, remainder};
}
/// @brief Gaussian-integer gcd, normalized to positive real part, or positive
/// imaginary part when real is zero; associates differ only by a unit.
inline gaussian_integer gaussian_gcd(gaussian_integer first,
gaussian_integer second) {
while (second != gaussian_integer{}) {
auto [quotient, remainder] = gaussian_divmod(first, second);
(void)quotient;
first = second;
second = remainder;
}
if (first.real < 0 || (first.real == 0 && first.imaginary < 0)) {
first = -first;
}
return first;
}
} // namespace noya