gaussian_integer.hpp¶
提供高斯整数的四则与欧几里得除法；适合 \(\mathbb Z[i]\) 上 gcd、二平方和等问题。
\[
\displaystyle \frac{a+bi}{c+di}
\]
Complexity: Time: \(O(\log \max (\lVert a\rVert , \lVert b\rVert))\) Euclidean steps. Space: \(O(1)\).
AC 记录：gcd_of_gaussian_integers。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(log max(norm(a),norm(b))) Euclidean steps.
/// Space: O(1).
#include <cassert>
#include <cstdint>
#include <limits>
#include <utility>
namespace noya {
/// @brief Signed 64-bit Gaussian integers. Arithmetic results, rounded
/// quotients, remainders and normalized gcd coordinates must fit int64_t.
/// Norms use unsigned 128-bit integers, including at the coordinate limits.
struct gaussian_integer {
std::int64_t re = 0;
std::int64_t im = 0;
friend bool operator==(const gaussian_integer &,
const gaussian_integer &) = default;
gaussian_integer &operator+=(const gaussian_integer &oth) {
re += oth.re;
im += oth.im;
return *this;
}
gaussian_integer &operator-=(const gaussian_integer &oth) {
re -= oth.re;
im -= oth.im;
return *this;
}
gaussian_integer &operator*=(const gaussian_integer &oth) {
__int128 nr = __int128(re) * oth.re - __int128(im) * oth.im;
__int128 nim = __int128(re) * oth.im + __int128(im) * oth.re;
re = std::int64_t(nr);
im = std::int64_t(nim);
return *this;
}
friend gaussian_integer operator+(gaussian_integer a,
const gaussian_integer &b) {
return a += b;
}
friend gaussian_integer operator-(gaussian_integer a,
const gaussian_integer &b) {
return a -= b;
}
friend gaussian_integer operator*(gaussian_integer a,
const gaussian_integer &b) {
return a *= b;
}
friend gaussian_integer operator-(gaussian_integer val) {
assert(val.re != std::numeric_limits<std::int64_t>::min() &&
val.im != std::numeric_limits<std::int64_t>::min());
return {-val.re, -val.im};
}
unsigned __int128 norm() const {
return static_cast<unsigned __int128>(__int128(re) * re) +
static_cast<unsigned __int128>(__int128(im) * im);
}
};
namespace gaussian_integer_internal {
using u128 = unsigned __int128;
inline std::int64_t round_div(u128 num, u128 den, bool neg) {
u128 q, r;
if ((num | den) >> 64 == 0) {
auto a = std::uint64_t(num), b = std::uint64_t(den);
q = a / b;
r = a % b;
} else {
q = num / den;
r = num % den;
}
if (r >= den - r) ++q; // Nearest integer; half ties go away from zero.
assert(q <= u128(std::numeric_limits<std::int64_t>::max()) + neg);
auto bits = std::uint64_t(q);
return std::int64_t(neg ? std::uint64_t(0) - bits : bits);
}
// Each term is a product of two int64_t coordinates. Their sum can be +2^127.
inline std::int64_t round_sum(__int128 a, __int128 b, u128 den) {
bool neg = a < 0;
u128 x = neg ? -u128(a) : u128(a);
u128 y = b < 0 ? -u128(b) : u128(b);
if (neg == (b < 0)) x += y;
else if (x >= y) x -= y;
else x = y - x, neg = !neg;
return round_div(x, den, neg);
}
} // namespace gaussian_integer_internal
inline std::int64_t gaussian_round_div(__int128 num, __int128 den) {
assert(den > 0);
using namespace gaussian_integer_internal;
return round_div(num < 0 ? -u128(num) : u128(num), u128(den), num < 0);
}
/// @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 dvd, gaussian_integer div) {
auto den = div.norm();
assert(den != 0);
using gaussian_integer_internal::round_sum;
gaussian_integer quo{
round_sum(__int128(dvd.re) * div.re, __int128(dvd.im) * div.im, den),
round_sum(__int128(dvd.im) * div.re, -__int128(dvd.re) * div.im, den)};
// The product may exceed int64_t even when both the quotient and remainder fit.
__int128 re = __int128(dvd.re) - __int128(quo.re) * div.re + __int128(quo.im) * div.im;
__int128 im = __int128(dvd.im) - __int128(quo.re) * div.im - __int128(quo.im) * div.re;
assert(re >= std::numeric_limits<std::int64_t>::min() &&
re <= std::numeric_limits<std::int64_t>::max() &&
im >= std::numeric_limits<std::int64_t>::min() &&
im <= std::numeric_limits<std::int64_t>::max());
gaussian_integer rem{std::int64_t(re), std::int64_t(im)};
assert(rem.norm() < den);
return {quo, rem};
}
/// @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 a, gaussian_integer b) {
while (b != gaussian_integer{}) {
auto [quo, rem] = gaussian_divmod(a, b);
(void)quo;
a = b;
b = rem;
}
if (a.re < 0 || (a.re == 0 && a.im < 0)) {
a = -a;
}
return a;
}
} // namespace noya
// NOYA_GAUSSIAN_INTEGER_HPP
#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 <limits>
#include <utility>
namespace noya {
/// @brief Signed 64-bit Gaussian integers. Arithmetic results, rounded
/// quotients, remainders and normalized gcd coordinates must fit int64_t.
/// Norms use unsigned 128-bit integers, including at the coordinate limits.
struct gaussian_integer {
std::int64_t re = 0;
std::int64_t im = 0;
friend bool operator==(const gaussian_integer &,
const gaussian_integer &) = default;
gaussian_integer &operator+=(const gaussian_integer &oth) {
re += oth.re;
im += oth.im;
return *this;
}
gaussian_integer &operator-=(const gaussian_integer &oth) {
re -= oth.re;
im -= oth.im;
return *this;
}
gaussian_integer &operator*=(const gaussian_integer &oth) {
__int128 nr = __int128(re) * oth.re - __int128(im) * oth.im;
__int128 nim = __int128(re) * oth.im + __int128(im) * oth.re;
re = std::int64_t(nr);
im = std::int64_t(nim);
return *this;
}
friend gaussian_integer operator+(gaussian_integer a,
const gaussian_integer &b) {
return a += b;
}
friend gaussian_integer operator-(gaussian_integer a,
const gaussian_integer &b) {
return a -= b;
}
friend gaussian_integer operator*(gaussian_integer a,
const gaussian_integer &b) {
return a *= b;
}
friend gaussian_integer operator-(gaussian_integer val) {
assert(val.re != std::numeric_limits<std::int64_t>::min() &&
val.im != std::numeric_limits<std::int64_t>::min());
return {-val.re, -val.im};
}
unsigned __int128 norm() const {
return static_cast<unsigned __int128>(__int128(re) * re) +
static_cast<unsigned __int128>(__int128(im) * im);
}
};
namespace gaussian_integer_internal {
using u128 = unsigned __int128;
inline std::int64_t round_div(u128 num, u128 den, bool neg) {
u128 q, r;
if ((num | den) >> 64 == 0) {
auto a = std::uint64_t(num), b = std::uint64_t(den);
q = a / b;
r = a % b;
} else {
q = num / den;
r = num % den;
}
if (r >= den - r) ++q; // Nearest integer; half ties go away from zero.
assert(q <= u128(std::numeric_limits<std::int64_t>::max()) + neg);
auto bits = std::uint64_t(q);
return std::int64_t(neg ? std::uint64_t(0) - bits : bits);
}
// Each term is a product of two int64_t coordinates. Their sum can be +2^127.
inline std::int64_t round_sum(__int128 a, __int128 b, u128 den) {
bool neg = a < 0;
u128 x = neg ? -u128(a) : u128(a);
u128 y = b < 0 ? -u128(b) : u128(b);
if (neg == (b < 0)) x += y;
else if (x >= y) x -= y;
else x = y - x, neg = !neg;
return round_div(x, den, neg);
}
} // namespace gaussian_integer_internal
inline std::int64_t gaussian_round_div(__int128 num, __int128 den) {
assert(den > 0);
using namespace gaussian_integer_internal;
return round_div(num < 0 ? -u128(num) : u128(num), u128(den), num < 0);
}
/// @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 dvd, gaussian_integer div) {
auto den = div.norm();
assert(den != 0);
using gaussian_integer_internal::round_sum;
gaussian_integer quo{
round_sum(__int128(dvd.re) * div.re, __int128(dvd.im) * div.im, den),
round_sum(__int128(dvd.im) * div.re, -__int128(dvd.re) * div.im, den)};
// The product may exceed int64_t even when both the quotient and remainder fit.
__int128 re = __int128(dvd.re) - __int128(quo.re) * div.re + __int128(quo.im) * div.im;
__int128 im = __int128(dvd.im) - __int128(quo.re) * div.im - __int128(quo.im) * div.re;
assert(re >= std::numeric_limits<std::int64_t>::min() &&
re <= std::numeric_limits<std::int64_t>::max() &&
im >= std::numeric_limits<std::int64_t>::min() &&
im <= std::numeric_limits<std::int64_t>::max());
gaussian_integer rem{std::int64_t(re), std::int64_t(im)};
assert(rem.norm() < den);
return {quo, rem};
}
/// @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 a, gaussian_integer b) {
while (b != gaussian_integer{}) {
auto [quo, rem] = gaussian_divmod(a, b);
(void)quo;
a = b;
b = rem;
}
if (a.re < 0 || (a.re == 0 && a.im < 0)) {
a = -a;
}
return a;
}
} // namespace noya
#endif // NOYA_GAUSSIAN_INTEGER_HPP
/// @complexity Time: O(log max(norm(a),norm(b))) Euclidean steps.
/// Space: O(1).
#include <cassert>
#include <cstdint>
#include <limits>
#include <utility>
namespace noya {
/// @brief Signed 64-bit Gaussian integers. Arithmetic results, rounded
/// quotients, remainders and normalized gcd coordinates must fit int64_t.
/// Norms use unsigned 128-bit integers, including at the coordinate limits.
struct gaussian_integer {
std::int64_t re = 0;
std::int64_t im = 0;
friend bool operator==(const gaussian_integer &,
const gaussian_integer &) = default;
gaussian_integer &operator+=(const gaussian_integer &oth) {
re += oth.re;
im += oth.im;
return *this;
}
gaussian_integer &operator-=(const gaussian_integer &oth) {
re -= oth.re;
im -= oth.im;
return *this;
}
gaussian_integer &operator*=(const gaussian_integer &oth) {
__int128 nr = __int128(re) * oth.re - __int128(im) * oth.im;
__int128 nim = __int128(re) * oth.im + __int128(im) * oth.re;
re = std::int64_t(nr);
im = std::int64_t(nim);
return *this;
}
friend gaussian_integer operator+(gaussian_integer a,
const gaussian_integer &b) {
return a += b;
}
friend gaussian_integer operator-(gaussian_integer a,
const gaussian_integer &b) {
return a -= b;
}
friend gaussian_integer operator*(gaussian_integer a,
const gaussian_integer &b) {
return a *= b;
}
friend gaussian_integer operator-(gaussian_integer val) {
assert(val.re != std::numeric_limits<std::int64_t>::min() &&
val.im != std::numeric_limits<std::int64_t>::min());
return {-val.re, -val.im};
}
unsigned __int128 norm() const {
return static_cast<unsigned __int128>(__int128(re) * re) +
static_cast<unsigned __int128>(__int128(im) * im);
}
};
namespace gaussian_integer_internal {
using u128 = unsigned __int128;
inline std::int64_t round_div(u128 num, u128 den, bool neg) {
u128 q, r;
if ((num | den) >> 64 == 0) {
auto a = std::uint64_t(num), b = std::uint64_t(den);
q = a / b;
r = a % b;
} else {
q = num / den;
r = num % den;
}
if (r >= den - r) ++q; // Nearest integer; half ties go away from zero.
assert(q <= u128(std::numeric_limits<std::int64_t>::max()) + neg);
auto bits = std::uint64_t(q);
return std::int64_t(neg ? std::uint64_t(0) - bits : bits);
}
// Each term is a product of two int64_t coordinates. Their sum can be +2^127.
inline std::int64_t round_sum(__int128 a, __int128 b, u128 den) {
bool neg = a < 0;
u128 x = neg ? -u128(a) : u128(a);
u128 y = b < 0 ? -u128(b) : u128(b);
if (neg == (b < 0)) x += y;
else if (x >= y) x -= y;
else x = y - x, neg = !neg;
return round_div(x, den, neg);
}
} // namespace gaussian_integer_internal
inline std::int64_t gaussian_round_div(__int128 num, __int128 den) {
assert(den > 0);
using namespace gaussian_integer_internal;
return round_div(num < 0 ? -u128(num) : u128(num), u128(den), num < 0);
}
/// @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 dvd, gaussian_integer div) {
auto den = div.norm();
assert(den != 0);
using gaussian_integer_internal::round_sum;
gaussian_integer quo{
round_sum(__int128(dvd.re) * div.re, __int128(dvd.im) * div.im, den),
round_sum(__int128(dvd.im) * div.re, -__int128(dvd.re) * div.im, den)};
// The product may exceed int64_t even when both the quotient and remainder fit.
__int128 re = __int128(dvd.re) - __int128(quo.re) * div.re + __int128(quo.im) * div.im;
__int128 im = __int128(dvd.im) - __int128(quo.re) * div.im - __int128(quo.im) * div.re;
assert(re >= std::numeric_limits<std::int64_t>::min() &&
re <= std::numeric_limits<std::int64_t>::max() &&
im >= std::numeric_limits<std::int64_t>::min() &&
im <= std::numeric_limits<std::int64_t>::max());
gaussian_integer rem{std::int64_t(re), std::int64_t(im)};
assert(rem.norm() < den);
return {quo, rem};
}
/// @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 a, gaussian_integer b) {
while (b != gaussian_integer{}) {
auto [quo, rem] = gaussian_divmod(a, b);
(void)quo;
a = b;
b = rem;
}
if (a.re < 0 || (a.re == 0 && a.im < 0)) {
a = -a;
}
return a;
}
} // namespace noya
// NOYA_GAUSSIAN_INTEGER_HPP