algo.5 Inverse of a real number (#15)

Author: tian-lt-personal

src/epsilon/chars.hpp

@@ -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 = coro::sync_get(num(n));
+    auto xn = num(n).get();
     auto sgn = xn.sgn;  // use the absolute value of xn to round towards zero.
     xn.sgn = sign::positive;
 

src/epsilon/coro.hpp

@@ -32,6 +32,8 @@ class lazy {
   class awaiter;
   awaiter operator co_await() noexcept;
 
+  T get();
+
  private:
   explicit lazy(std::coroutine_handle<lazy_promise<T>> coro) : coro_(coro) {}
   std::coroutine_handle<lazy_promise<T>> coro_;
@@ -46,7 +48,14 @@ class lazy_promise {
   constexpr auto initial_suspend() noexcept { return std::suspend_always{}; }
   constexpr auto final_suspend() noexcept {
     struct final_awaiter : std::suspend_always {
-      auto await_suspend(std::coroutine_handle<lazy_promise> handle) noexcept { return handle.promise().cont_; }
+      std::coroutine_handle<> await_suspend(std::coroutine_handle<lazy_promise> handle) noexcept {
+        auto cont = handle.promise().cont_;
+        if (cont) {
+          return cont;
+        } else {
+          return std::noop_coroutine();
+        }
+      }
     };
     return final_awaiter{};
   }
@@ -73,13 +82,18 @@ class lazy<T>::awaiter : public std::suspend_always {
   T await_resume() {
     return std::visit(
         tmp::overloads{
-            [](T& val) -> T&& { return std::move(val); },                                 // normal path
-            [](std::monostate) -> T&& { std::terminate(); },                              // unassigned
-            [](std::exception_ptr& ex) -> T&& { std::rethrow_exception(std::move(ex)); }  // contains exception
+            [](T& val) -> T { return std::move(val); },                                 // normal path
+            [](std::monostate) -> T { std::terminate(); },                              // unassigned
+            [](std::exception_ptr& ex) -> T { std::rethrow_exception(std::move(ex)); }  // contains exception
         },
         coro_.promise().state_);
   }
 
+  T get() {
+    coro_();
+    return await_resume();
+  }
+
  private:
   std::coroutine_handle<lazy_promise<T>> coro_;
 };
@@ -88,15 +102,13 @@ typename lazy<T>::awaiter lazy<T>::operator co_await() noexcept {
   return lazy<T>::awaiter{coro_};
 }
 
-struct forget {};
-
 template <class T>
-T sync_get(lazy<T> awaitable) {
-  std::optional<T> slot;
-  [&]() -> forget { slot = co_await awaitable; }();
-  return std::move(*slot);
+T lazy<T>::get() {
+  return lazy<T>::awaiter{coro_}.get();
 }
 
+struct forget {};
+
 }  // namespace epx::coro
 
 namespace std {

src/epsilon/r.hpp

@@ -7,6 +7,7 @@
 // std
 #include <algorithm>
 #include <functional>
+#include <limits>
 
 // epx
 #include "coro.hpp"
@@ -73,11 +74,20 @@ coro::lazy<int> msd(r<C> x, int max) {
       if (c > 0 || i >= max) {
         co_return i;
       }
-      xi = co_await x(++i);  // TODO: use faster search algo.
+      ++i;  // TODO: use faster search algo.
+      if (i > max_msd<global_config_tag>) [[unlikely]] {
+        throw msd_overflow_error{};
+      }
+      xi = co_await x(i);
     }
   }
 }
 
+template <container C>
+coro::lazy<int> msd(r<C> x) {
+  return msd(std::move(x), std::numeric_limits<int>::max());
+}
+
 template <container C>
 constexpr r<C> mul(r<C> x, r<C> y) {
   struct {
@@ -100,6 +110,28 @@ constexpr r<C> mul(r<C> x, r<C> y) {
   return expr;
 }
 
+template <container C>
+constexpr r<C> inv(r<C> x) {
+  struct {
+    r<C> x;
+
+    coro::lazy<z<C>> operator()(unsigned int n) {
+      int msdx = co_await msd(x);
+      int nn = static_cast<int>(n);
+      if (nn <= -msdx) {
+        co_return z<C>{};  // zero
+      }
+
+      int k = nn + 2 * msdx + 1;
+      auto p = mul_4exp(details::one<C>(), k + nn);
+      auto q = add(co_await x(k), details::one<C>());
+      auto [quo, _] = ceil_div(p, q);
+      co_return add(quo, details::one<C>());
+    }
+  } expr{.x = std::move(x)};
+  return expr;
+}
+
 }  // namespace epx
 
 #endif  // EPSILON_INC_R_HPP

src/epsilon/t.hpp

@@ -9,16 +9,16 @@
 
 namespace epx {
 
-using default_digit_t = uint32_t;
-using max_digit_t = uint32_t;
-using default_container_t = std::vector<default_digit_t>;
+using default_digit_type = uint32_t;
+using max_digit_type = uint32_t;
+using default_container_type = std::vector<default_digit_type>;
 
-struct bad_digit_t;
+struct bad_digit_type;
 template <class D>
-using wide_digit_t =
+using wide_digit_type =
     std::conditional_t<sizeof(D) == sizeof(uint8_t), uint16_t,
                        std::conditional_t<sizeof(D) == sizeof(uint16_t), uint32_t,
-                                          std::conditional_t<sizeof(D) == sizeof(uint32_t), uint64_t, bad_digit_t>>>;
+                                          std::conditional_t<sizeof(D) == sizeof(uint32_t), uint64_t, bad_digit_type>>>;
 
 template <class T>
 concept container = std::ranges::random_access_range<T> &&  //
@@ -29,12 +29,12 @@ concept container = std::ranges::random_access_range<T> &&  //
                     } && requires(T c) {
                       c.push_back(typename T::value_type{});
                       c.reserve(size_t{});
-                      sizeof(typename T::value_type) < sizeof(max_digit_t);
+                      sizeof(typename T::value_type) < sizeof(max_digit_type);
                     };
 
 enum class sign : uint8_t { positive, negative };
 
-template <container C = default_container_t>
+template <container C = default_container_type>
 struct z {
   using container_type = C;
   using digit_type = typename C::value_type;
@@ -44,6 +44,11 @@ struct z {
   sign sgn = sign::positive;
 };
 
+struct global_config_tag {};  // allow user-defined configs
+
+template <typename>
+constexpr int max_msd = 10000;  // can be overridden by global_config_tag
+
 struct divide_by_zero_error : public std::runtime_error {
   divide_by_zero_error() : std::runtime_error("epx: divide by zero") {}
 };
@@ -52,6 +57,10 @@ struct overflow_error : public std::runtime_error {
   overflow_error() : std::runtime_error("epx: overflow error") {}
 };
 
+struct msd_overflow_error : public std::runtime_error {
+  msd_overflow_error() : std::runtime_error("epx: msd overflow error") {}
+};
+
 }  // namespace epx
 
 #endif

src/epsilon/z.hpp

@@ -76,7 +76,7 @@ constexpr auto umul(T lhs, T rhs) {
 // u0 - LSD, u1 - MSD
 template <class D>
 constexpr auto div_2d(D u0, D u1, D v) {
-  using W = wide_digit_t<D>;
+  using W = wide_digit_type<D>;
   struct result_t {
     W q;
     D r;
@@ -233,7 +233,7 @@ constexpr z<C> mul_n(const z<C>& lhs, const z<C>& rhs) {
   for (size_t j = 0; j < std::ranges::size(b); ++j) {
     D cy = 0;
     for (size_t i = 0; i < std::ranges::size(a); ++i) {
-      auto [p0, p1] = details::umul<default_digit_t>(a[i], b[j]);
+      auto [p0, p1] = details::umul<default_digit_type>(a[i], b[j]);
       p0 += cy;
       cy = (p0 < cy ? 1u : 0u) + p1;
       r.digits[i + j] += p0;
@@ -248,7 +248,7 @@ constexpr z<C> mul_n(const z<C>& lhs, const z<C>& rhs) {
 template <container C>
 constexpr auto div_n(const z<C>& u, typename z<C>::digit_type v) {
   using D = typename z<C>::digit_type;
-  static_assert(sizeof(D) <= sizeof(max_digit_t));
+  static_assert(sizeof(D) <= sizeof(max_digit_type));
 
   struct {
     z<C> q;
@@ -283,7 +283,7 @@ constexpr auto div_n(z<C> lhs, z<C> rhs) {
   if (rel > 0) {
     if (std::ranges::size(rhs.digits) > 1) {
       using D = typename z<C>::digit_type;
-      using W = wide_digit_t<D>;
+      using W = wide_digit_type<D>;
       constexpr W b = W{1} << (sizeof(D) * CHAR_BIT);
       constexpr D mask = D{b - 1};
 
@@ -426,6 +426,26 @@ constexpr auto floor_div(z<C> lhs, z<C> rhs) {
   return res;
 }
 
+template <container C>
+constexpr auto ceil_div(z<C> lhs, z<C> rhs) {
+  struct result_t {
+    z<C> q;
+    z<C> r;
+  };
+
+  auto sgn = lhs.sgn == rhs.sgn ? sign::positive : sign::negative;
+  auto [q, r] = div_n(std::move(lhs), rhs);
+  result_t res = {.q = std::move(q), .r = std::move(r)};
+
+  if (sgn == sign::positive && !is_zero(res.r)) {
+    res.q = add_n(res.q, details::one<C>());
+    res.r = sub_n(rhs, res.r);
+  }
+  res.q.sgn = sgn;
+  res.r.sgn = rhs.sgn == sign::positive ? sign::negative : sign::positive;
+  return res;
+}
+
 template <container C>
 constexpr z<C>& mul_2exp(z<C>& val, int exp) {
   using D = typename z<C>::digit_type;