gaussian_integer.hpp¶
提供高斯整数的四则与欧几里得除法；适合 \(\mathbb Z[i]\) 上 gcd、二平方和等问题。
\[
\displaystyle \\frac{a+bi}{c+di}
\]
Complexity: Time: O(log max(norm(a),norm(b))) 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 <cstdlib>
#include <utility>
namespace noya {
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) {
return {-val.re, -val.im};
}
__int128 norm() const { return __int128(re) * re + __int128(im) * im; }
};
inline std::int64_t gaussian_round_div(__int128 num, __int128 den) {
assert(den > 0);
__int128 quo = num / den;
__int128 rem = num % den;
__int128 ar = rem < 0 ? -rem : rem;
if (2 * ar >= den) {
quo += num >= 0 ? 1 : -1;
}
return std::int64_t(quo);
}
/// @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) {
__int128 den = div.norm();
assert(den != 0);
__int128 rn = __int128(dvd.re) * div.re + __int128(dvd.im) * div.im;
__int128 ni = __int128(dvd.im) * div.re - __int128(dvd.re) * div.im;
gaussian_integer quo{gaussian_round_div(rn, den),
gaussian_round_div(ni, den)};
gaussian_integer rem = dvd - quo * div;
assert(rem.norm() < div.norm());
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
#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 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) {
return {-val.re, -val.im};
}
__int128 norm() const { return __int128(re) * re + __int128(im) * im; }
};
inline std::int64_t gaussian_round_div(__int128 num, __int128 den) {
assert(den > 0);
__int128 quo = num / den;
__int128 rem = num % den;
__int128 ar = rem < 0 ? -rem : rem;
if (2 * ar >= den) {
quo += num >= 0 ? 1 : -1;
}
return std::int64_t(quo);
}
/// @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) {
__int128 den = div.norm();
assert(den != 0);
__int128 rn = __int128(dvd.re) * div.re + __int128(dvd.im) * div.im;
__int128 ni = __int128(dvd.im) * div.re - __int128(dvd.re) * div.im;
gaussian_integer quo{gaussian_round_div(rn, den),
gaussian_round_div(ni, den)};
gaussian_integer rem = dvd - quo * div;
assert(rem.norm() < div.norm());
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
#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 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) {
return {-val.re, -val.im};
}
__int128 norm() const { return __int128(re) * re + __int128(im) * im; }
};
inline std::int64_t gaussian_round_div(__int128 num, __int128 den) {
assert(den > 0);
__int128 quo = num / den;
__int128 rem = num % den;
__int128 ar = rem < 0 ? -rem : rem;
if (2 * ar >= den) {
quo += num >= 0 ? 1 : -1;
}
return std::int64_t(quo);
}
/// @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) {
__int128 den = div.norm();
assert(den != 0);
__int128 rn = __int128(dvd.re) * div.re + __int128(dvd.im) * div.im;
__int128 ni = __int128(dvd.im) * div.re - __int128(dvd.re) * div.im;
gaussian_integer quo{gaussian_round_div(rn, den),
gaussian_round_div(ni, den)};
gaussian_integer rem = dvd - quo * div;
assert(rem.norm() < div.norm());
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