rational_reconstruction.hpp¶
Reconstruct a reduced fraction numerator/denominator congruent to a residue modulo modulus within supplied bounds. Uniqueness is guaranteed when 2numerator_bounddenominator_bound < modulus.
\[
\displaystyle \frac{x}{y}\equiv r\pmod m
\]
Implementation¶
#ifndef NOYA_RATIONAL_RECONSTRUCTION_HPP
#define NOYA_RATIONAL_RECONSTRUCTION_HPP 1
/// @complexity Time: O(log m).
/// Space: O(1).
#include "noya/fraction.hpp"
#include <cassert>
#include <cstdlib>
#include <cstdint>
#include <numeric>
#include <optional>
namespace noya {
/// @brief Reconstruct a reduced fraction numerator/denominator congruent to a
/// residue modulo modulus within supplied bounds. Uniqueness is guaranteed when
/// 2*numerator_bound*denominator_bound < modulus.
inline std::optional<fraction>
rational_reconstruction(std::int64_t residue, std::int64_t modulus,
std::int64_t numerator_bound,
std::int64_t denominator_bound) {
assert(modulus >= 2);
assert(numerator_bound >= 0 && denominator_bound >= 1);
residue %= modulus;
if (residue < 0) {
residue += modulus;
}
std::int64_t old_remainder = modulus;
std::int64_t remainder = residue;
std::int64_t old_denominator = 0;
std::int64_t denominator = 1;
while (std::llabs(remainder) > numerator_bound) {
if (remainder == 0) {
return std::nullopt;
}
std::int64_t quotient = old_remainder / remainder;
__int128 next_remainder =
__int128(old_remainder) - __int128(quotient) * remainder;
__int128 next_denominator =
__int128(old_denominator) - __int128(quotient) * denominator;
old_remainder = remainder;
remainder = std::int64_t(next_remainder);
old_denominator = denominator;
denominator = std::int64_t(next_denominator);
}
if (denominator < 0) {
denominator = -denominator;
remainder = -remainder;
}
if (denominator == 0 || denominator > denominator_bound ||
std::gcd(remainder, denominator) != 1 ||
std::gcd(denominator, modulus) != 1) {
return std::nullopt;
}
__int128 congruence = __int128(residue) * denominator - remainder;
if (congruence % modulus != 0) {
return std::nullopt;
}
return fraction(remainder, denominator);
}
} // namespace noya
#endif // NOYA_RATIONAL_RECONSTRUCTION_HPP
#include <cassert>
#include <compare>
#include <cstdint>
#include <cstdlib>
#include <numeric>
#include <optional>
/// @complexity Time: O(log m).
/// Space: O(1).
/// @complexity Time: O(log max(|numerator|,denominator)) normalization; arithmetic itself is O(1).
/// Space: O(1).
namespace noya {
/// @brief Normalized signed 64-bit rational number; arithmetic results must fit
/// int64, while comparison and intermediate products use signed 128-bit.
struct fraction {
std::int64_t numerator = 0;
std::int64_t denominator = 1;
fraction() = default;
fraction(std::int64_t numerator_, std::int64_t denominator_ = 1)
: numerator(numerator_), denominator(denominator_) {
normalize();
}
void normalize() {
assert(denominator != 0);
if (denominator < 0) {
numerator = -numerator;
denominator = -denominator;
}
std::int64_t divisor = std::gcd(numerator, denominator);
numerator /= divisor;
denominator /= divisor;
}
friend bool operator==(const fraction &, const fraction &) = default;
friend std::strong_ordering operator<=>(const fraction &first,
const fraction &second) {
__int128 left = __int128(first.numerator) * second.denominator;
__int128 right = __int128(second.numerator) * first.denominator;
return left < right ? std::strong_ordering::less
: left > right ? std::strong_ordering::greater
: std::strong_ordering::equal;
}
friend fraction operator+(const fraction &first, const fraction &second) {
std::int64_t divisor = std::gcd(first.denominator, second.denominator);
__int128 numerator =
__int128(first.numerator) * (second.denominator / divisor) +
__int128(second.numerator) * (first.denominator / divisor);
__int128 denominator =
__int128(first.denominator / divisor) * second.denominator;
return {std::int64_t(numerator), std::int64_t(denominator)};
}
friend fraction operator-(const fraction &first, const fraction &second) {
return first + fraction(-second.numerator, second.denominator);
}
friend fraction operator*(fraction first, fraction second) {
std::int64_t first_divisor =
std::gcd(first.numerator < 0 ? -first.numerator : first.numerator,
second.denominator);
std::int64_t second_divisor =
std::gcd(second.numerator < 0 ? -second.numerator : second.numerator,
first.denominator);
first.numerator /= first_divisor;
second.denominator /= first_divisor;
second.numerator /= second_divisor;
first.denominator /= second_divisor;
return {std::int64_t(__int128(first.numerator) * second.numerator),
std::int64_t(__int128(first.denominator) * second.denominator)};
}
friend fraction operator/(const fraction &first, const fraction &second) {
assert(second.numerator != 0);
return first * fraction(second.denominator, second.numerator);
}
friend fraction operator-(const fraction &value) {
return {-value.numerator, value.denominator};
}
};
} // namespace noya
namespace noya {
/// @brief Reconstruct a reduced fraction numerator/denominator congruent to a
/// residue modulo modulus within supplied bounds. Uniqueness is guaranteed when
/// 2*numerator_bound*denominator_bound < modulus.
inline std::optional<fraction>
rational_reconstruction(std::int64_t residue, std::int64_t modulus,
std::int64_t numerator_bound,
std::int64_t denominator_bound) {
assert(modulus >= 2);
assert(numerator_bound >= 0 && denominator_bound >= 1);
residue %= modulus;
if (residue < 0) {
residue += modulus;
}
std::int64_t old_remainder = modulus;
std::int64_t remainder = residue;
std::int64_t old_denominator = 0;
std::int64_t denominator = 1;
while (std::llabs(remainder) > numerator_bound) {
if (remainder == 0) {
return std::nullopt;
}
std::int64_t quotient = old_remainder / remainder;
__int128 next_remainder =
__int128(old_remainder) - __int128(quotient) * remainder;
__int128 next_denominator =
__int128(old_denominator) - __int128(quotient) * denominator;
old_remainder = remainder;
remainder = std::int64_t(next_remainder);
old_denominator = denominator;
denominator = std::int64_t(next_denominator);
}
if (denominator < 0) {
denominator = -denominator;
remainder = -remainder;
}
if (denominator == 0 || denominator > denominator_bound ||
std::gcd(remainder, denominator) != 1 ||
std::gcd(denominator, modulus) != 1) {
return std::nullopt;
}
__int128 congruence = __int128(residue) * denominator - remainder;
if (congruence % modulus != 0) {
return std::nullopt;
}
return fraction(remainder, denominator);
}
} // namespace noya