rational_reconstruction.hpp¶
从模意义下的剩余值，在给定分子分母界内唯一恢复最简分数。
\[
\displaystyle \frac{x}{y}\equiv r\pmod m
\]
Complexity: Time: O(log m). Space: O(1).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(log m).
/// Space: O(1).
#include "noya/math/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 res modulo mod within supplied bounds. Uniqueness is guaranteed when
/// 2*nb*db < mod.
inline std::optional<fraction> rational_reconstruction(std::int64_t res,
std::int64_t mod,
std::int64_t nb,
std::int64_t db) {
assert(mod >= 2);
assert(nb >= 0 && db >= 1);
res %= mod;
if (res < 0) {
res += mod;
}
std::int64_t r0 = mod;
std::int64_t rem = res;
std::int64_t d0 = 0;
std::int64_t den = 1;
while (std::llabs(rem) > nb) {
if (rem == 0) {
return std::nullopt;
}
std::int64_t quo = r0 / rem;
__int128 r1 = __int128(r0) - __int128(quo) * rem;
__int128 d1 = __int128(d0) - __int128(quo) * den;
r0 = rem;
rem = std::int64_t(r1);
d0 = den;
den = std::int64_t(d1);
}
if (den < 0) {
den = -den;
rem = -rem;
}
if (den == 0 || den > db || std::gcd(rem, den) != 1 ||
std::gcd(den, mod) != 1) {
return std::nullopt;
}
__int128 con = __int128(res) * den - rem;
if (con % mod != 0) {
return std::nullopt;
}
return fraction(rem, den);
}
} // namespace noya
// NOYA_RATIONAL_RECONSTRUCTION_HPP
#ifndef NOYA_RATIONAL_RECONSTRUCTION_HPP
#define NOYA_RATIONAL_RECONSTRUCTION_HPP 1
/// @complexity Time: O(log m).
/// Space: O(1).
#include "noya/math/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 res modulo mod within supplied bounds. Uniqueness is guaranteed when
/// 2*nb*db < mod.
inline std::optional<fraction> rational_reconstruction(std::int64_t res,
std::int64_t mod,
std::int64_t nb,
std::int64_t db) {
assert(mod >= 2);
assert(nb >= 0 && db >= 1);
res %= mod;
if (res < 0) {
res += mod;
}
std::int64_t r0 = mod;
std::int64_t rem = res;
std::int64_t d0 = 0;
std::int64_t den = 1;
while (std::llabs(rem) > nb) {
if (rem == 0) {
return std::nullopt;
}
std::int64_t quo = r0 / rem;
__int128 r1 = __int128(r0) - __int128(quo) * rem;
__int128 d1 = __int128(d0) - __int128(quo) * den;
r0 = rem;
rem = std::int64_t(r1);
d0 = den;
den = std::int64_t(d1);
}
if (den < 0) {
den = -den;
rem = -rem;
}
if (den == 0 || den > db || std::gcd(rem, den) != 1 ||
std::gcd(den, mod) != 1) {
return std::nullopt;
}
__int128 con = __int128(res) * den - rem;
if (con % mod != 0) {
return std::nullopt;
}
return fraction(rem, den);
}
} // namespace noya
#endif // NOYA_RATIONAL_RECONSTRUCTION_HPP
/// @complexity Time: O(log m).
/// Space: O(1).
/// @complexity Time: O(log max(|num|,den)) normalization; arithmetic itself is O(1).
/// Space: O(1).
#include <cassert>
#include <compare>
#include <cstdint>
#include <limits>
#include <numeric>
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 num = 0;
std::int64_t den = 1;
fraction() = default;
fraction(std::int64_t nx, std::int64_t dx = 1) : num(nx), den(dx) {
normalize();
}
void normalize() {
*this = from_wide(num, den);
}
friend bool operator==(const fraction &, const fraction &) = default;
friend std::strong_ordering operator<=>(const fraction &a,
const fraction &b) {
__int128 l = __int128(a.num) * b.den;
__int128 r = __int128(b.num) * a.den;
return l < r ? std::strong_ordering::less
: l > r ? std::strong_ordering::greater
: std::strong_ordering::equal;
}
friend fraction operator+(const fraction &a, const fraction &b) {
std::int64_t div = std::gcd(a.den, b.den);
__int128 num =
__int128(a.num) * (b.den / div) + __int128(b.num) * (a.den / div);
__int128 den = __int128(a.den / div) * b.den;
return from_wide(num, den);
}
friend fraction operator-(const fraction &a, const fraction &b) {
return from_wide(__int128(a.num) * b.den - __int128(b.num) * a.den,
__int128(a.den) * b.den);
}
friend fraction operator*(fraction a, fraction b) {
return from_wide(__int128(a.num) * b.num, __int128(a.den) * b.den);
}
friend fraction operator/(const fraction &a, const fraction &b) {
assert(b.num != 0);
return from_wide(__int128(a.num) * b.den, __int128(a.den) * b.num);
}
friend fraction operator-(const fraction &val) {
return from_wide(-__int128(val.num), val.den);
}
private:
// Reduce before narrowing: an intermediate numerator or denominator may
// exceed int64 even when the normalized result is representable.
static fraction from_wide(__int128 n, __int128 d) {
assert(d != 0);
if (d < 0) {
n = -n;
d = -d;
}
__int128 a = n < 0 ? -n : n, b = d;
while (b != 0) {
__int128 rem = a % b;
a = b;
b = rem;
}
n /= a;
d /= a;
assert(n >= std::numeric_limits<std::int64_t>::lowest() &&
n <= std::numeric_limits<std::int64_t>::max() &&
d <= std::numeric_limits<std::int64_t>::max());
fraction res;
res.num = std::int64_t(n);
res.den = std::int64_t(d);
return res;
}
};
} // namespace noya
#include <cassert>
#include <cstdlib>
#include <cstdint>
#include <numeric>
#include <optional>
namespace noya {
/// @brief Reconstruct a reduced fraction numerator/denominator congruent to a
/// residue res modulo mod within supplied bounds. Uniqueness is guaranteed when
/// 2*nb*db < mod.
inline std::optional<fraction> rational_reconstruction(std::int64_t res,
std::int64_t mod,
std::int64_t nb,
std::int64_t db) {
assert(mod >= 2);
assert(nb >= 0 && db >= 1);
res %= mod;
if (res < 0) {
res += mod;
}
std::int64_t r0 = mod;
std::int64_t rem = res;
std::int64_t d0 = 0;
std::int64_t den = 1;
while (std::llabs(rem) > nb) {
if (rem == 0) {
return std::nullopt;
}
std::int64_t quo = r0 / rem;
__int128 r1 = __int128(r0) - __int128(quo) * rem;
__int128 d1 = __int128(d0) - __int128(quo) * den;
r0 = rem;
rem = std::int64_t(r1);
d0 = den;
den = std::int64_t(d1);
}
if (den < 0) {
den = -den;
rem = -rem;
}
if (den == 0 || den > db || std::gcd(rem, den) != 1 ||
std::gcd(den, mod) != 1) {
return std::nullopt;
}
__int128 con = __int128(res) * den - rem;
if (con % mod != 0) {
return std::nullopt;
}
return fraction(rem, den);
}
} // namespace noya
// NOYA_RATIONAL_RECONSTRUCTION_HPP