Skip to content

refactor the real number type r #17

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tian-lt-personal merged 2 commits into from
Feb 18, 2026
Merged

refactor the real number type r #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7ca6683
mpa
tian-lt-personal Feb 18, 2026
d347516
update test case
tian-lt-personal Feb 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/epsilon/chars.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ constexpr std::string to_string(r<C> num, unsigned int k) {
constexpr int extra_precision = 10;
if constexpr (B == 10) {
auto n = static_cast<int>(log_4_10 * k) + extra_precision;
auto xn = num(n).get();
auto xn = num.approx(n).get();
auto sgn = xn.sgn; // use the absolute value of xn to round towards zero.
xn.sgn = sign::positive;

Expand Down
113 changes: 60 additions & 53 deletions src/epsilon/r.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,37 +16,55 @@
namespace epx {

template <container C>
using r = std::function<coro::lazy<z<C>>(int)>; // q (n) = floor(p * 4^n / q)
class r {
public:
template <class F>
explicit r(F&& x) : x_(std::forward<F>(x)) {}
coro::lazy<z<C>> approx(int n) const {
if (n <= mpa_) {
co_return mul_4exp(const_cast<const z<C>&>(x_mpa_), n - mpa_);
} else {
mpa_ = n;
x_mpa_ = co_await x_(n);
co_return x_mpa_;
}
}

private:
std::function<coro::lazy<z<C>>(int)> x_;
mutable int mpa_ = std::numeric_limits<int>::min(); // most precise approximation
mutable z<C> x_mpa_;
};

template <container C>
constexpr r<C> make_q(z<C> p, z<C> q) {
return [p = std::move(p), q = std::move(q)](int n) -> coro::lazy<z<C>> {
return r<C>{[p = std::move(p), q = std::move(q)](int n) -> coro::lazy<z<C>> {
auto rr = mul_4exp(p, n);
auto [quo, _] = floor_div(rr, q);
co_return quo;
};
co_return quo; // q (n) = floor(p * 4^n / q)
}};
}

template <container C>
constexpr r<C> add(r<C> x, r<C> y) {
return [x = std::move(x), y = std::move(y)](int n) -> coro::lazy<z<C>> {
return r<C>{[x = std::move(x), y = std::move(y)](int n) -> coro::lazy<z<C>> {
if (n == std::numeric_limits<int>::max()) [[unlikely]] {
throw precision_overflow_error{};
}
// w = 4 when B = 4.
auto xn = co_await x(n + 1);
auto yn = co_await y(n + 1);
auto xn = co_await x.approx(n + 1);
auto yn = co_await y.approx(n + 1);
co_return mul_4exp(xn + yn, -1);
};
}};
}

template <container C>
constexpr r<C> opp(r<C> x) {
return [x = std::move(x)](int n) -> coro::lazy<z<C>> {
auto xn = co_await x(n);
return r<C>{[x = std::move(x)](int n) -> coro::lazy<z<C>> {
auto xn = co_await x.approx(n);
negate(xn);
co_return xn;
};
}};
}

template <container C>
Expand All @@ -55,11 +73,11 @@ coro::lazy<int> msd(r<C> x, int max) {
const auto one = details::one<C>();
const auto base = details::base<C>();

auto x0 = co_await x(0);
auto x0 = co_await x.approx(0);
if (cmp_n(x0, base) > 0) { // 4 < x0
int i = -1;
for (;;) {
auto xi = co_await x(i);
auto xi = co_await x.approx(i);
if (cmp_n(xi, one) <= 0) {
co_return i + 1;
}
Expand All @@ -81,7 +99,7 @@ coro::lazy<int> msd(r<C> x, int max) {
if (i > max_msd<global_config_tag>) [[unlikely]] {
throw msd_overflow_error{};
}
xi = co_await x(i);
xi = co_await x.approx(i);
}
}
}
Expand All @@ -93,52 +111,41 @@ coro::lazy<int> msd(r<C> x) {

template <container C>
constexpr r<C> mul(r<C> x, r<C> y) {
struct {
r<C> x;
r<C> y;

coro::lazy<z<C>> operator()(int n) {
if (n > std::numeric_limits<int>::max() - 3) [[unlikely]] {
throw precision_overflow_error{};
}
const auto one = details::one<C>();
int max_msd = n + 3 - (n + 2) / 2;
int px = n - (co_await msd(y, max_msd)) + 3;
int py = n - (co_await msd(x, max_msd)) + 3;
auto xpx = co_await x(px);
auto ypy = co_await y(py);
auto xy = mul_4exp(add_n(mul_n(xpx, ypy), one), n - px - py);
xy.sgn = xpx.sgn == ypy.sgn ? sign::positive : sign::negative;
co_return xy;
return r<C>{[x = std::move(x), y = std::move(y)](int n) -> coro::lazy<z<C>> {
if (n > std::numeric_limits<int>::max() - 3) [[unlikely]] {
throw precision_overflow_error{};
}
} expr{.x = std::move(x), .y = std::move(y)};
return expr;
const auto one = details::one<C>();
int max_msd = n + 3 - (n + 2) / 2;
int px = n - (co_await msd(y, max_msd)) + 3;
int py = n - (co_await msd(x, max_msd)) + 3;
auto xpx = co_await x.approx(px);
auto ypy = co_await y.approx(py);
auto xy = mul_4exp(add_n(mul_n(xpx, ypy), one), n - px - py);
xy.sgn = xpx.sgn == ypy.sgn ? sign::positive : sign::negative;
co_return xy;
}};
}

template <container C>
constexpr r<C> inv(r<C> x) {
struct {
r<C> x;

coro::lazy<z<C>> operator()(int n) {
int msdx = co_await msd(x);
if (n <= -msdx) {
co_return z<C>{}; // zero
}

assert(std::numeric_limits<int>::max() > 2 * max_msd<global_config_tag> + 1);
if (n > std::numeric_limits<int>::max() - 2 * max_msd<global_config_tag> - 1) [[unlikely]] {
throw precision_overflow_error{};
}
return r<C>{[x = std::move(x)](int n) -> coro::lazy<z<C>> {
int msdx = co_await msd(x);
if (n <= -msdx) {
co_return z<C>{}; // zero
}

int k = n + 2 * msdx + 1;
auto p = mul_4exp(details::one<C>(), k + n);
auto q = add(co_await x(k), details::one<C>());
auto [quo, _] = ceil_div(p, q);
co_return add(quo, details::one<C>());
assert(std::numeric_limits<int>::max() > 2 * max_msd<global_config_tag> + 1);
if (n > std::numeric_limits<int>::max() - 2 * max_msd<global_config_tag> - 1) [[unlikely]] {
throw precision_overflow_error{};
}
} expr{.x = std::move(x)};
return expr;

int k = n + 2 * msdx + 1;
auto p = mul_4exp(details::one<C>(), k + n);
auto q = add(co_await x.approx(k), details::one<C>());
auto [quo, _] = ceil_div(p, q);
co_return add(quo, details::one<C>());
}};
}

} // namespace epx
Expand Down
22 changes: 11 additions & 11 deletions src/ut/r_tests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,28 +18,28 @@ TEST(r_tests, make_q) {
{
sz zero;
auto q = epx::make_q(zero, one); // 0/1
EXPECT_TRUE(epx::is_zero(q(1).get()));
EXPECT_TRUE(epx::is_zero(q(10).get()));
EXPECT_TRUE(epx::is_zero(q(100).get()));
EXPECT_TRUE(epx::is_zero(q.approx(1).get()));
EXPECT_TRUE(epx::is_zero(q.approx(10).get()));
EXPECT_TRUE(epx::is_zero(q.approx(100).get()));
}
{
auto q = epx::make_q(one, one); // 1/1
EXPECT_EQ(stosz("4"), q(1).get());
EXPECT_EQ(stosz("16"), q(2).get());
EXPECT_EQ(stosz("1048576"), q(10).get());
EXPECT_EQ(stosz("4"), q.approx(1).get());
EXPECT_EQ(stosz("16"), q.approx(2).get());
EXPECT_EQ(stosz("1048576"), q.approx(10).get());
}
{
auto num = stosz("3");
auto q = epx::make_q(one, num); // 1/3
EXPECT_EQ(one, q(1).get());
EXPECT_EQ(stosz("341"), q(5).get());
EXPECT_EQ(one, q.approx(1).get());
EXPECT_EQ(stosz("341"), q.approx(5).get());
}

{
auto num = stosz("-3");
auto q = epx::make_q(one, num); // 1/-3
EXPECT_EQ(stosz("-2"), q(1).get());
EXPECT_EQ(stosz("-342"), q(5).get());
EXPECT_EQ(stosz("-2"), q.approx(1).get());
EXPECT_EQ(stosz("-342"), q.approx(5).get());
}
}

Expand Down Expand Up @@ -210,7 +210,7 @@ TEST(r_tests, inv) {
{
auto zero = epx::make_q(stosz("0"), stosz("1"));
auto expr = epx::inv(zero);
EXPECT_THROW(expr(10).get(), epx::msd_overflow_error);
EXPECT_THROW(expr.approx(10).get(), epx::msd_overflow_error);
}
{
auto r = epx::make_q(stosz("1"), stosz("1"));
Expand Down