rational_reconstruction.hpp¶
从模剩余值恢复最简分数；分子绝对值上界、分母上界和模数分别为 \(N,D,m\)。当 \(2ND<m\) 时保证找到唯一解或判定无解；超出此范围可能漏解。
\[
\displaystyle x\equiv ry\pmod m,\quad |x|\le N,\quad 1\le y\le D,\quad\gcd(y,m)=1
\]
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 x/y with |x| <= nb, 1 <= y <= db, gcd(y,mod) == 1,
/// and x == res*y (mod mod). If 2*nb*db < mod, returns the unique reduced
/// fraction when it exists, otherwise nullopt. With larger bounds, a returned
/// fraction is still valid, but nullopt does not prove that no solution exists.
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 x/y with |x| <= nb, 1 <= y <= db, gcd(y,mod) == 1,
/// and x == res*y (mod mod). If 2*nb*db < mod, returns the unique reduced
/// fraction when it exists, otherwise nullopt. With larger bounds, a returned
/// fraction is still valid, but nullopt does not prove that no solution exists.
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)) construction/arithmetic; O(1) comparison.
/// 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() {
assert(den != 0);
auto g = std::gcd(magnitude(num), magnitude(den));
*this = from_reduced(__int128(num) / g, __int128(den) / g);
}
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) {
return add(a, b, 1);
}
friend fraction operator-(const fraction &a, const fraction &b) {
return add(a, b, -1);
}
friend fraction operator*(const fraction &a, const fraction &b) {
auto g = std::gcd(magnitude(a.num), std::uint64_t(b.den));
auto h = std::gcd(magnitude(b.num), std::uint64_t(a.den));
return from_reduced(__int128(a.num / std::int64_t(g)) * (b.num / std::int64_t(h)),
__int128(a.den / std::int64_t(h)) * (b.den / std::int64_t(g)));
}
friend fraction operator/(const fraction &a, const fraction &b) {
assert(b.num != 0);
auto g = std::gcd(magnitude(a.num), magnitude(b.num));
auto h = std::gcd(a.den, b.den);
return from_reduced((__int128(a.num) / g) * (b.den / h),
(__int128(b.num) / g) * (a.den / h));
}
friend fraction operator-(const fraction &val) {
return from_reduced(-__int128(val.num), val.den);
}
private:
static std::uint64_t magnitude(std::int64_t n) {
return n < 0 ? std::uint64_t(0) - std::uint64_t(n) : std::uint64_t(n);
}
static fraction add(const fraction &a, const fraction &b, int sgn) {
auto g = std::gcd(a.den, b.den);
auto ad = a.den / g, bd = b.den / g;
__int128 n = __int128(a.num) * bd + sgn * __int128(b.num) * ad;
// Coprime ad and bd leave only factors of g to cancel.
auto h = std::gcd(std::int64_t(n % g), g);
return from_reduced(n / h, __int128(ad) * (b.den / h));
}
static fraction from_reduced(__int128 n, __int128 d) {
assert(d != 0);
if (d < 0) {
n = -n;
d = -d;
}
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 x/y with |x| <= nb, 1 <= y <= db, gcd(y,mod) == 1,
/// and x == res*y (mod mod). If 2*nb*db < mod, returns the unique reduced
/// fraction when it exists, otherwise nullopt. With larger bounds, a returned
/// fraction is still valid, but nullopt does not prove that no solution exists.
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