Sum of Two Squares
(math/two_square_sum.hpp)
- View this file on GitHub
- Last update: 2026-07-15 01:09:58+09:00
- Include:
#include "math/two_square_sum.hpp"
Overview
Find representations of an integer as a sum of two non-negative squares.
#include "math/two_square_sum.hpp"
All names are in m1une::math.
For example, 25 has four ordered non-negative representations:
0^2 + 5^2
3^2 + 4^2
4^2 + 3^2
5^2 + 0^2
so two_square_sum(25) returns the pairs (0, 5), (3, 4), (4, 3), and
(5, 0).
API
std::vector<std::pair<uint64_t, uint64_t>> two_square_sum(uint64_t value);
std::vector<std::pair<uint64_t, uint64_t>>
represent_as_two_square_sum(uint64_t value);
bool is_two_square_sum(uint64_t value);
| Function | Description |
|---|---|
two_square_sum(value) |
Returns all ordered non-negative pairs (a, b) such that value = a * a + b * b. |
represent_as_two_square_sum(value) |
Alias of two_square_sum(value). |
is_two_square_sum(value) |
Returns whether at least one representation exists. |
The returned vector is sorted in lexicographic order and contains no duplicate
pairs. Both (a, b) and (b, a) are included when they are different.
two_square_sum(0) returns one pair: (0, 0).
Criterion
By Fermat’s theorem on sums of two squares, an integer has a representation as
$a^2 + b^2$ exactly when every prime congruent to 3 modulo 4 appears with an
even exponent in its prime factorization.
This is the condition tested by is_two_square_sum.
For construction, each prime congruent to 1 modulo 4 is first represented as
a Gaussian integer norm, and all choices of conjugate factors are combined.
Prime powers congruent to 3 modulo 4 only contribute to a real scale factor
because their exponents must be even.
Complexity
The function first factors value with Pollard-Rho from
prime_factorization.hpp. After factorization, the construction is linear in
the number of generated Gaussian-factor choices and the number of returned
pairs.
For a prime factorization
\[value = 2^e \prod p_i^{\alpha_i} \prod q_j^{2\beta_j},\]where every p_i is 1 modulo 4 and every q_j is 3 modulo 4, the
number of Gaussian-factor choices is
The returned vector contains at most twice this many pairs before duplicate removal.
Example
#include "math/two_square_sum.hpp"
#include <cstdint>
#include <iostream>
#include <utility>
#include <vector>
int main() {
uint64_t value = 65;
std::vector<std::pair<uint64_t, uint64_t>> pairs =
m1une::math::two_square_sum(value);
for (const auto& pair : pairs) {
std::cout << pair.first << " " << pair.second << "\n";
}
}
Depends on
- Integer Roots and Powers
(math/integer_arithmetic.hpp)
- 64-bit Prime Factorization
(math/prime_factorization.hpp)
Required by
Verified with
Code
#ifndef M1UNE_MATH_TWO_SQUARE_SUM_HPP
#define M1UNE_MATH_TWO_SQUARE_SUM_HPP 1
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <utility>
#include <vector>
#include "integer_arithmetic.hpp"
#include "prime_factorization.hpp"
namespace m1une {
namespace math {
namespace two_square_sum_detail {
struct GaussianInteger {
__int128_t real;
__int128_t imaginary;
};
inline GaussianInteger multiply(GaussianInteger first, GaussianInteger second) {
return GaussianInteger{first.real * second.real - first.imaginary * second.imaginary,
first.real * second.imaginary + first.imaginary * second.real};
}
inline GaussianInteger power(GaussianInteger base, int exponent) {
GaussianInteger result{1, 0};
while (exponent > 0) {
if ((exponent & 1) != 0) result = multiply(result, base);
exponent >>= 1;
if (exponent != 0) base = multiply(base, base);
}
return result;
}
inline uint64_t absolute_value(__int128_t value) {
return static_cast<uint64_t>(value < 0 ? -value : value);
}
inline uint64_t pow_uint64(uint64_t base, int exponent) {
uint64_t result = 1;
while (exponent > 0) {
if ((exponent & 1) != 0) result *= base;
exponent >>= 1;
if (exponent != 0) base *= base;
}
return result;
}
inline GaussianInteger prime_one_mod_four_representation(uint64_t prime) {
assert(prime % 4 == 1);
uint64_t non_residue = 2;
while (internal::power_mod(non_residue, (prime - 1) / 2, prime) != prime - 1) {
non_residue++;
}
uint64_t root = internal::power_mod(non_residue, (prime - 1) / 4, prime);
uint64_t previous = prime;
uint64_t current = root;
while (static_cast<__uint128_t>(current) * current > prime) {
uint64_t next = previous % current;
previous = current;
current = next;
}
const uint64_t real = current;
const uint64_t remaining = prime - real * real;
const uint64_t imaginary = isqrt(remaining);
assert(imaginary * imaginary == remaining);
return GaussianInteger{static_cast<__int128_t>(real), static_cast<__int128_t>(imaginary)};
}
inline std::vector<GaussianInteger> prime_power_choices(GaussianInteger factor, int exponent) {
GaussianInteger conjugate{factor.real, -factor.imaginary};
std::vector<GaussianInteger> positive_powers(exponent + 1);
std::vector<GaussianInteger> negative_powers(exponent + 1);
positive_powers[0] = GaussianInteger{1, 0};
negative_powers[0] = GaussianInteger{1, 0};
for (int i = 0; i < exponent; i++) {
positive_powers[i + 1] = multiply(positive_powers[i], factor);
negative_powers[i + 1] = multiply(negative_powers[i], conjugate);
}
std::vector<GaussianInteger> result;
result.reserve(exponent + 1);
for (int take_positive = 0; take_positive <= exponent; take_positive++) {
result.push_back(multiply(positive_powers[take_positive], negative_powers[exponent - take_positive]));
}
return result;
}
inline void add_nonnegative_pairs(std::vector<std::pair<uint64_t, uint64_t>>& result, GaussianInteger value,
uint64_t scale) {
const uint64_t first = absolute_value(value.real) * scale;
const uint64_t second = absolute_value(value.imaginary) * scale;
result.emplace_back(first, second);
result.emplace_back(second, first);
}
} // namespace two_square_sum_detail
inline std::vector<std::pair<uint64_t, uint64_t>> two_square_sum(uint64_t value) {
using two_square_sum_detail::GaussianInteger;
if (value == 0) return std::vector<std::pair<uint64_t, uint64_t>>{std::pair<uint64_t, uint64_t>(0, 0)};
uint64_t real_scale = 1;
std::vector<std::vector<GaussianInteger>> choices_by_prime;
for (const auto& factor : prime_factorize(value)) {
const uint64_t prime = factor.first;
const int exponent = factor.second;
if (prime % 4 == 3) {
if ((exponent & 1) != 0) return {};
real_scale *= two_square_sum_detail::pow_uint64(prime, exponent / 2);
continue;
}
GaussianInteger gaussian_factor;
if (prime == 2) {
gaussian_factor = GaussianInteger{1, 1};
} else {
gaussian_factor = two_square_sum_detail::prime_one_mod_four_representation(prime);
}
choices_by_prime.push_back(two_square_sum_detail::prime_power_choices(gaussian_factor, exponent));
}
std::vector<GaussianInteger> values;
values.push_back(GaussianInteger{1, 0});
for (const std::vector<GaussianInteger>& choices : choices_by_prime) {
std::vector<GaussianInteger> next;
next.reserve(values.size() * choices.size());
for (GaussianInteger value_so_far : values) {
for (GaussianInteger choice : choices) {
next.push_back(two_square_sum_detail::multiply(value_so_far, choice));
}
}
values.swap(next);
}
std::vector<std::pair<uint64_t, uint64_t>> result;
result.reserve(values.size() * 2);
for (GaussianInteger gaussian : values) {
two_square_sum_detail::add_nonnegative_pairs(result, gaussian, real_scale);
}
std::sort(result.begin(), result.end());
result.erase(std::unique(result.begin(), result.end()), result.end());
return result;
}
inline std::vector<std::pair<uint64_t, uint64_t>> represent_as_two_square_sum(uint64_t value) {
return two_square_sum(value);
}
inline bool is_two_square_sum(uint64_t value) {
if (value == 0) return true;
for (const auto& factor : prime_factorize(value)) {
if (factor.first % 4 == 3 && (factor.second & 1) != 0) return false;
}
return true;
}
} // namespace math
} // namespace m1une
#endif // M1UNE_MATH_TWO_SQUARE_SUM_HPP#line 1 "math/two_square_sum.hpp"
#include <algorithm>
#include <cassert>
#include <cstdint>
#include <utility>
#include <vector>
#line 1 "math/integer_arithmetic.hpp"
#line 5 "math/integer_arithmetic.hpp"
#include <concepts>
#include <limits>
#include <optional>
#include <type_traits>
namespace m1une {
namespace math {
namespace integer_arithmetic_detail {
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr std::optional<T> checked_multiply(T first, T second) {
constexpr T minimum = std::numeric_limits<T>::min();
constexpr T maximum = std::numeric_limits<T>::max();
if constexpr (std::unsigned_integral<T>) {
if (second != 0 && maximum / second < first) return std::nullopt;
} else {
if (0 < first) {
if (0 < second) {
if (maximum / second < first) return std::nullopt;
} else if (second < minimum / first) {
return std::nullopt;
}
} else if (first < 0) {
if (0 < second) {
if (first < minimum / second) return std::nullopt;
} else if (second < maximum / first) {
return std::nullopt;
}
}
}
return T(first * second);
}
template <std::unsigned_integral T>
constexpr bool kth_power_leq(T base, unsigned exponent, T limit) {
assert(exponent > 0);
if (base <= 1) return base <= limit;
const T multiplication_limit = limit / base;
T product = 1;
for (unsigned i = 0; i < exponent; i++) {
if (product > multiplication_limit) return false;
product *= base;
}
return true;
}
} // namespace integer_arithmetic_detail
// Returns floor(sqrt(value)) exactly, without floating-point arithmetic.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T isqrt(T value) {
if constexpr (std::signed_integral<T>) assert(0 <= value);
if (value <= 1) return value;
T low = 1;
T high = value / 2 + 1;
while (low < high) {
T middle = low + (high - low + 1) / 2;
if (middle <= value / middle) {
low = middle;
} else {
high = middle - 1;
}
}
return low;
}
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T floor_sqrt(T value) {
return isqrt(value);
}
// Returns ceil(sqrt(value)) exactly, without floating-point arithmetic.
template <std::integral T>
requires(!std::same_as<std::remove_cv_t<T>, bool>)
constexpr T ceil_sqrt(T value) {
T result = isqrt(value);
if (result == 0) return 0;
if (result != 0 && value / result == result && value % result == 0) {
return result;
}
return result + 1;
}
// Returns floor(value^(1 / degree)) exactly, without floating-point arithmetic.
template <std::integral T, std::integral Degree>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Degree>, bool>
)
constexpr T floor_kth_root(T value, Degree degree) {
if constexpr (std::signed_integral<T>) {
assert(0 <= value);
if (value < 0) return T();
}
assert(0 < degree);
if (degree <= 0) return T();
if (value <= 1 || degree == 1) return value;
if (degree == 2) return isqrt(value);
using U = std::make_unsigned_t<T>;
using UDegree = std::make_unsigned_t<Degree>;
constexpr int digits = std::numeric_limits<U>::digits;
const UDegree unsigned_degree = static_cast<UDegree>(degree);
if (unsigned_degree >= static_cast<UDegree>(digits)) return T(1);
const unsigned exponent = static_cast<unsigned>(unsigned_degree);
const U unsigned_value = static_cast<U>(value);
int bit_width = 0;
for (U remaining = unsigned_value; remaining != 0; remaining >>= 1) {
bit_width++;
}
const int root_bits =
(bit_width + static_cast<int>(exponent) - 1) /
static_cast<int>(exponent);
U low = 1;
U high = U(1) << root_bits;
while (high - low > 1) {
const U middle = low + (high - low) / 2;
if (
integer_arithmetic_detail::kth_power_leq(
middle, exponent, unsigned_value
)
) {
low = middle;
} else {
high = middle;
}
}
return static_cast<T>(low);
}
// Returns base^exponent, or nullopt when the result does not fit in T.
template <std::integral T, std::unsigned_integral Exponent>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr std::optional<T> checked_ipow(T base, Exponent exponent) {
T result = 1;
while (exponent != 0) {
if (exponent & 1) {
auto product =
integer_arithmetic_detail::checked_multiply(result, base);
if (!product.has_value()) return std::nullopt;
result = *product;
}
exponent >>= 1;
if (exponent != 0) {
auto square =
integer_arithmetic_detail::checked_multiply(base, base);
if (!square.has_value()) return std::nullopt;
base = *square;
}
}
return result;
}
template <std::integral T, std::unsigned_integral Exponent>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr std::optional<T> checked_integer_pow(T base, Exponent exponent) {
return checked_ipow(base, exponent);
}
// Returns base^exponent. The result must be representable by T.
template <std::integral T, std::unsigned_integral Exponent>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr T ipow(T base, Exponent exponent) {
std::optional<T> result = checked_ipow(base, exponent);
assert(result.has_value());
return result.value_or(T());
}
template <std::integral T, std::unsigned_integral Exponent>
requires(
!std::same_as<std::remove_cv_t<T>, bool>
&& !std::same_as<std::remove_cv_t<Exponent>, bool>
)
constexpr T integer_pow(T base, Exponent exponent) {
return ipow(base, exponent);
}
} // namespace math
} // namespace m1une
#line 1 "math/prime_factorization.hpp"
#line 7 "math/prime_factorization.hpp"
#include <numeric>
#line 10 "math/prime_factorization.hpp"
namespace m1une {
namespace math {
namespace internal {
inline uint64_t multiply_mod(uint64_t a, uint64_t b, uint64_t mod) {
return static_cast<uint64_t>(static_cast<unsigned __int128>(a) * b % mod);
}
inline uint64_t power_mod(uint64_t base, uint64_t exponent, uint64_t mod) {
uint64_t result = 1;
while (exponent > 0) {
if (exponent & 1) result = multiply_mod(result, base, mod);
base = multiply_mod(base, base, mod);
exponent >>= 1;
}
return result;
}
inline uint64_t pollard_random() {
static uint64_t state = 0x123456789abcdef0ULL;
state += 0x9e3779b97f4a7c15ULL;
uint64_t value = state;
value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL;
value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL;
return value ^ (value >> 31);
}
} // namespace internal
inline bool is_prime(uint64_t value) {
if (value < 2) return false;
for (uint64_t prime : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
if (value % prime == 0) return value == prime;
}
uint64_t odd_part = value - 1;
int power_of_two = 0;
while ((odd_part & 1) == 0) {
odd_part >>= 1;
power_of_two++;
}
for (uint64_t base : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
if (base % value == 0) continue;
uint64_t x = internal::power_mod(base % value, odd_part, value);
if (x == 1 || x == value - 1) continue;
bool composite = true;
for (int i = 1; i < power_of_two; i++) {
x = internal::multiply_mod(x, x, value);
if (x == value - 1) {
composite = false;
break;
}
}
if (composite) return false;
}
return true;
}
namespace internal {
inline uint64_t pollard_rho(uint64_t value) {
for (uint64_t prime : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
if (value % prime == 0) return prime;
}
while (true) {
const uint64_t constant = pollard_random() % (value - 1) + 1;
uint64_t y = pollard_random() % (value - 1) + 1;
uint64_t x = 0;
uint64_t saved_y = 0;
uint64_t gcd = 1;
uint64_t segment_length = 1;
auto advance = [&](uint64_t current) {
return static_cast<uint64_t>(
(static_cast<unsigned __int128>(multiply_mod(current, current, value)) + constant) % value);
};
while (gcd == 1) {
x = y;
for (uint64_t i = 0; i < segment_length; i++) y = advance(y);
for (uint64_t offset = 0; offset < segment_length && gcd == 1; offset += 128) {
saved_y = y;
uint64_t product = 1;
const uint64_t block = std::min<uint64_t>(128, segment_length - offset);
for (uint64_t i = 0; i < block; i++) {
y = advance(y);
const uint64_t difference = x > y ? x - y : y - x;
product = multiply_mod(product, difference, value);
}
gcd = std::gcd(product, value);
}
segment_length <<= 1;
}
if (gcd == value) {
do {
saved_y = advance(saved_y);
const uint64_t difference = x > saved_y ? x - saved_y : saved_y - x;
gcd = std::gcd(difference, value);
} while (gcd == 1);
}
if (gcd != value) return gcd;
}
}
inline void factor_recursively(uint64_t value, std::vector<uint64_t>& factors) {
if (value == 1) return;
if (is_prime(value)) {
factors.push_back(value);
return;
}
const uint64_t divisor = pollard_rho(value);
factor_recursively(divisor, factors);
factor_recursively(value / divisor, factors);
}
} // namespace internal
inline std::vector<uint64_t> prime_factors(uint64_t value) {
assert(value >= 1);
std::vector<uint64_t> result;
internal::factor_recursively(value, result);
std::sort(result.begin(), result.end());
return result;
}
inline std::vector<std::pair<uint64_t, int>> prime_factorize(uint64_t value) {
std::vector<uint64_t> factors = prime_factors(value);
std::vector<std::pair<uint64_t, int>> result;
for (uint64_t prime : factors) {
if (result.empty() || result.back().first != prime) {
result.emplace_back(prime, 1);
} else {
result.back().second++;
}
}
return result;
}
inline std::vector<uint64_t> divisors(uint64_t value) {
std::vector<uint64_t> result = {1};
for (const auto& factor : prime_factorize(value)) {
const int current_size = int(result.size());
uint64_t power = 1;
for (int exponent = 1; exponent <= factor.second; exponent++) {
power *= factor.first;
for (int i = 0; i < current_size; i++) {
result.push_back(result[i] * power);
}
}
}
std::sort(result.begin(), result.end());
return result;
}
inline uint64_t euler_phi(uint64_t value) {
assert(value >= 1);
uint64_t result = value;
for (const auto& factor : prime_factorize(value)) {
result = result / factor.first * (factor.first - 1);
}
return result;
}
inline int mobius(uint64_t value) {
assert(value >= 1);
int result = 1;
for (const auto& factor : prime_factorize(value)) {
if (factor.second >= 2) return 0;
result = -result;
}
return result;
}
} // namespace math
} // namespace m1une
#line 12 "math/two_square_sum.hpp"
namespace m1une {
namespace math {
namespace two_square_sum_detail {
struct GaussianInteger {
__int128_t real;
__int128_t imaginary;
};
inline GaussianInteger multiply(GaussianInteger first, GaussianInteger second) {
return GaussianInteger{first.real * second.real - first.imaginary * second.imaginary,
first.real * second.imaginary + first.imaginary * second.real};
}
inline GaussianInteger power(GaussianInteger base, int exponent) {
GaussianInteger result{1, 0};
while (exponent > 0) {
if ((exponent & 1) != 0) result = multiply(result, base);
exponent >>= 1;
if (exponent != 0) base = multiply(base, base);
}
return result;
}
inline uint64_t absolute_value(__int128_t value) {
return static_cast<uint64_t>(value < 0 ? -value : value);
}
inline uint64_t pow_uint64(uint64_t base, int exponent) {
uint64_t result = 1;
while (exponent > 0) {
if ((exponent & 1) != 0) result *= base;
exponent >>= 1;
if (exponent != 0) base *= base;
}
return result;
}
inline GaussianInteger prime_one_mod_four_representation(uint64_t prime) {
assert(prime % 4 == 1);
uint64_t non_residue = 2;
while (internal::power_mod(non_residue, (prime - 1) / 2, prime) != prime - 1) {
non_residue++;
}
uint64_t root = internal::power_mod(non_residue, (prime - 1) / 4, prime);
uint64_t previous = prime;
uint64_t current = root;
while (static_cast<__uint128_t>(current) * current > prime) {
uint64_t next = previous % current;
previous = current;
current = next;
}
const uint64_t real = current;
const uint64_t remaining = prime - real * real;
const uint64_t imaginary = isqrt(remaining);
assert(imaginary * imaginary == remaining);
return GaussianInteger{static_cast<__int128_t>(real), static_cast<__int128_t>(imaginary)};
}
inline std::vector<GaussianInteger> prime_power_choices(GaussianInteger factor, int exponent) {
GaussianInteger conjugate{factor.real, -factor.imaginary};
std::vector<GaussianInteger> positive_powers(exponent + 1);
std::vector<GaussianInteger> negative_powers(exponent + 1);
positive_powers[0] = GaussianInteger{1, 0};
negative_powers[0] = GaussianInteger{1, 0};
for (int i = 0; i < exponent; i++) {
positive_powers[i + 1] = multiply(positive_powers[i], factor);
negative_powers[i + 1] = multiply(negative_powers[i], conjugate);
}
std::vector<GaussianInteger> result;
result.reserve(exponent + 1);
for (int take_positive = 0; take_positive <= exponent; take_positive++) {
result.push_back(multiply(positive_powers[take_positive], negative_powers[exponent - take_positive]));
}
return result;
}
inline void add_nonnegative_pairs(std::vector<std::pair<uint64_t, uint64_t>>& result, GaussianInteger value,
uint64_t scale) {
const uint64_t first = absolute_value(value.real) * scale;
const uint64_t second = absolute_value(value.imaginary) * scale;
result.emplace_back(first, second);
result.emplace_back(second, first);
}
} // namespace two_square_sum_detail
inline std::vector<std::pair<uint64_t, uint64_t>> two_square_sum(uint64_t value) {
using two_square_sum_detail::GaussianInteger;
if (value == 0) return std::vector<std::pair<uint64_t, uint64_t>>{std::pair<uint64_t, uint64_t>(0, 0)};
uint64_t real_scale = 1;
std::vector<std::vector<GaussianInteger>> choices_by_prime;
for (const auto& factor : prime_factorize(value)) {
const uint64_t prime = factor.first;
const int exponent = factor.second;
if (prime % 4 == 3) {
if ((exponent & 1) != 0) return {};
real_scale *= two_square_sum_detail::pow_uint64(prime, exponent / 2);
continue;
}
GaussianInteger gaussian_factor;
if (prime == 2) {
gaussian_factor = GaussianInteger{1, 1};
} else {
gaussian_factor = two_square_sum_detail::prime_one_mod_four_representation(prime);
}
choices_by_prime.push_back(two_square_sum_detail::prime_power_choices(gaussian_factor, exponent));
}
std::vector<GaussianInteger> values;
values.push_back(GaussianInteger{1, 0});
for (const std::vector<GaussianInteger>& choices : choices_by_prime) {
std::vector<GaussianInteger> next;
next.reserve(values.size() * choices.size());
for (GaussianInteger value_so_far : values) {
for (GaussianInteger choice : choices) {
next.push_back(two_square_sum_detail::multiply(value_so_far, choice));
}
}
values.swap(next);
}
std::vector<std::pair<uint64_t, uint64_t>> result;
result.reserve(values.size() * 2);
for (GaussianInteger gaussian : values) {
two_square_sum_detail::add_nonnegative_pairs(result, gaussian, real_scale);
}
std::sort(result.begin(), result.end());
result.erase(std::unique(result.begin(), result.end()), result.end());
return result;
}
inline std::vector<std::pair<uint64_t, uint64_t>> represent_as_two_square_sum(uint64_t value) {
return two_square_sum(value);
}
inline bool is_two_square_sum(uint64_t value) {
if (value == 0) return true;
for (const auto& factor : prime_factorize(value)) {
if (factor.first % 4 == 3 && (factor.second & 1) != 0) return false;
}
return true;
}
} // namespace math
} // namespace m1une