algo.4 multiplication of two real numbers

src/epsilon/r.hpp

@@ -5,6 +5,7 @@
 #define EPSILON_INC_R_HPP
 
 // std
+#include <algorithm>
 #include <functional>
 
 // epx
@@ -44,6 +45,61 @@ constexpr r<C> opp(r<C> x) {
   };
 }
 
+template <container C>
+coro::lazy<int> msd(r<C> x, int max) {
+  const auto zero = details::zero<C>();
+  const auto one = details::one<C>();
+  const auto base = details::base<C>();
+
+  auto x0 = co_await x(0);
+  if (cmp_n(x0, base) > 0) {  // 4 < x0
+    int i = -1;
+    for (;;) {
+      auto xi = co_await x(i);
+      if (cmp_n(xi, one) <= 0) {
+        co_return i + 1;
+      }
+      --i;  // TODO: use faster search algo.
+    }
+  } else if (cmp_n(x0, one) > 0) {  // 1 < x0 <= 4
+    co_return 0;
+  } else if (cmp_n(x0, zero) > 0) {  // 0 < x0 <= 1
+    co_return 1;
+  } else {  // x0 <= 0
+    int i = 0;
+    auto xi = x0;
+    for (;;) {
+      auto c = cmp_n(xi, one);
+      if (c > 0 || i >= max) {
+        co_return i;
+      }
+      xi = co_await x(++i);  // TODO: use faster search algo.
+    }
+  }
+}
+
+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()(unsigned int n) {
+      const auto one = details::one<C>();
+      int nn = static_cast<int>(n);
+      int max_msd = nn + 3 - (nn + 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;
+    }
+  } expr{.x = std::move(x), .y = std::move(y)};
+  return expr;
+}
+
 }  // namespace epx
 
 #endif  // EPSILON_INC_R_HPP

src/epsilon/z.hpp

@@ -29,6 +29,11 @@ constexpr z<C> one() {
   return z<C>{.digits = {1}};
 }
 
+template <container C>
+constexpr z<C> base() {
+  return z<C>{.digits = {4}};
+}
+
 template <container C>
 constexpr z<C> ten() {
   return z<C>{.digits = {10}};

src/ut/r_tests.cpp

@@ -118,4 +118,112 @@ TEST(r_tests, opp) {
   }
 }
 
+TEST(r_tests, msd) {
+  {
+    auto zero = epx::make_q(stosz("0"), stosz("1"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(10, msd);
+    msd = epx::coro::sync_get(epx::msd(zero, 1000));
+    EXPECT_EQ(1000, msd);
+  }
+  {
+    auto zero = epx::make_q(stosz("1"), stosz("1"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(1, msd);
+  }
+  {
+    auto zero = epx::make_q(stosz("11"), stosz("10"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(1, msd);
+  }
+  {
+    auto zero = epx::make_q(stosz("2"), stosz("1"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(0, msd);
+  }
+  {
+    auto zero = epx::make_q(stosz("1"), stosz("2"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(1, msd);
+  }
+  {
+    auto zero = epx::make_q(stosz("5"), stosz("1"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(0, msd);
+  }
+  {
+    auto zero = epx::make_q(stosz("10"), stosz("1"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(-1, msd);
+  }
+  {
+    auto zero = epx::make_q(stosz("127"), stosz("1"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(-2, msd);
+  }
+  {
+    auto zero = epx::make_q(stosz("128"), stosz("1"));
+    int msd = epx::coro::sync_get(epx::msd(zero, 10));
+    EXPECT_EQ(-3, msd);
+  }
+}
+
+TEST(r_tests, mul) {
+  {
+    auto x = epx::make_q(stosz("0"), stosz("1"));
+    auto y = epx::make_q(stosz("0"), stosz("1"));
+    auto expr = epx::mul(x, y);
+    EXPECT_EQ(
+        "0."
+        "00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"
+        "00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"
+        "00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"
+        "00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"
+        "000000000000000000000000000000000000000000000000000000000000",
+        epx::to_string(expr, 500));
+  }
+  {
+    auto x = epx::make_q(stosz("1"), stosz("1"));
+    auto y = epx::make_q(stosz("1"), stosz("1"));
+    auto expr = epx::mul(x, y);
+    EXPECT_EQ("1.00000000000000000000", epx::to_string(expr, 20));
+  }
+  {
+    auto x = epx::make_q(stosz("-2"), stosz("1"));
+    auto y = epx::make_q(stosz("1"), stosz("1"));
+    auto expr = epx::mul(x, y);
+    EXPECT_EQ("-2.0000000000000000000000000000000000000000", epx::to_string(expr, 40));
+  }
+  {
+    auto x = epx::make_q(stosz("1"), stosz("3"));
+    auto y = epx::make_q(stosz("3"), stosz("1"));
+    auto expr = epx::mul(x, y);
+    EXPECT_EQ("1.0000000000000000000000000000000000000000", epx::to_string(expr, 40));
+  }
+  {
+    auto x = epx::make_q(stosz("7"), stosz("5"));
+    auto y = epx::make_q(stosz("2"), stosz("9"));
+    auto expr = epx::mul(x, y);
+    EXPECT_EQ("0.3111111111111111111111111111111111111111", epx::to_string(expr, 40));
+  }
+  {
+    auto x = epx::make_q(stosz("11"), stosz("7"));
+    auto y = epx::make_q(stosz("1"), stosz("121"));
+    auto expr = epx::mul(x, y);
+    EXPECT_EQ("0.0129870129870129870129870129870129870130", epx::to_string(expr, 40));
+  }
+  {
+    auto x = epx::make_q(stosz("1"), stosz("10000000000"));
+    auto y = epx::make_q(stosz("1"), stosz("2"));
+    auto expr = epx::mul(x, y);
+    EXPECT_EQ("0.00000000005000000000", epx::to_string(expr, 20));
+  }
+  {
+    auto x = epx::make_q(stosz("827368917649287346"), stosz("1"));
+    auto y = epx::make_q(stosz("1"), stosz("792873649187263413"));
+    auto expr = epx::mul(x, y);
+    EXPECT_EQ("1.04350664005214875176", epx::to_string(expr, 20));
+  }
+}
+
 }  // namespace epxut