algo.6 kth root.

src/epsilon/chars.hpp

@@ -107,9 +107,9 @@ constexpr std::string to_string(r<C> num, unsigned int k) {
     d.sgn = sgn;      // restore the sign
 
     auto s = to_string(d);
-    size_t len = is_positive(d) ? s.length() : s.length() - 1;
+    size_t len = is_negative(d) ? s.length() - 1 : s.length();
     if (len <= k) {
-      s.insert(is_positive(d) ? 0 : 1, k - len + 1, '0');
+      s.insert(is_negative(d) ? 1 : 0, k - len + 1, '0');
     }
     if (k > 0) {
       s.insert(s.length() - k, 1, '.');

src/epsilon/r.hpp

@@ -148,6 +148,21 @@ constexpr r<C> inv(r<C> x) {
   }};
 }
 
+template <container C>
+constexpr r<C> root(r<C> x, int k) {
+  return r<C>{[x = std::move(x), k](int n) -> coro::lazy<z<C>> {
+    if (k < 2) [[unlikely]] {
+      throw kthroot_too_small_error{};
+    }
+
+    auto xn = co_await x.approx(n * k);
+    if (is_negative(xn)) [[unlikely]] {
+      throw negative_radicand_error{};
+    }
+    co_return root(xn, k);
+  }};
+}
+
 }  // namespace epx
 
 #endif  // EPSILON_INC_R_HPP

src/epsilon/t.hpp

@@ -61,6 +61,18 @@ struct precision_overflow_error : public std::runtime_error {
   precision_overflow_error() : std::runtime_error("epx: precision overflow error") {}
 };
 
+struct negative_radicand_error : public std::runtime_error {
+  negative_radicand_error() : std::runtime_error("epx: negative radicand error") {}
+};
+
+struct negative_zpow_error : public std::runtime_error {
+  negative_zpow_error() : std::runtime_error("epx: negative power for z error") {}
+};
+
+struct kthroot_too_small_error : public std::runtime_error {
+  kthroot_too_small_error() : std::runtime_error("epx: kth-root too small error") {}
+};
+
 }  // namespace epx
 
 #endif

src/epsilon/z.hpp

@@ -117,14 +117,36 @@ constexpr auto bit_shift(C& digits, int offset) {
 }
 }  // namespace details
 
+template <container C, std::integral T>
+constexpr z<C> create(T val) {
+  z<C> res;
+  if (val == 0) {
+    return res;
+  }
+  if constexpr (std::is_signed_v<T>) {
+    res.sgn = val < 0 ? sign::negative : sign::positive;
+    val = val < 0 ? -val : val;
+  }
+  using D = typename z<C>::digit_type;
+  if constexpr (sizeof(T) <= sizeof(D)) {
+    res.digits.push_back(static_cast<D>(val));
+  } else {
+    while (val > 0) {
+      res.digits.push_back(static_cast<D>(val));
+      val >>= (sizeof(D) * CHAR_BIT);
+    }
+  }
+  return res;
+}
+
 template <container C>
 constexpr bool is_zero(const z<C>& num) noexcept {
   return std::ranges::size(num.digits) == 0;
 }
 
 template <container C>
-constexpr bool is_positive(const z<C>& num) noexcept {
-  return num.sgn == sign::positive;
+constexpr bool is_negative(const z<C>& num) noexcept {
+  return num.sgn == sign::negative;
 }
 
 template <container C>
@@ -141,7 +163,7 @@ constexpr z<C>& negate(z<C>& num) noexcept {
   if (is_zero(num)) {
     return num;
   }
-  num.sgn = is_positive(num) ? sign::negative : sign::positive;
+  num.sgn = is_negative(num) ? sign::positive : sign::negative;
   return num;
 }
 
@@ -503,6 +525,50 @@ constexpr z<C> mul_4exp(const z<C>& val, int exp) {
   return mul_4exp(num, exp);
 }
 
+template <container C>
+constexpr z<C> pow(const z<C>& num, int exp) {
+  assert(exp >= 0);
+  z<C> res = details::one<C>();
+  if (exp > 0) {
+    for (int i = 0; i < exp; ++i) {
+      res = mul_n(res, num);
+    }
+  } else if (exp == 0) {
+    return res;
+  } else [[unlikely]] {
+    throw negative_zpow_error{};
+  }
+  if (is_negative(num) && (exp % 2 == 1)) {
+    res.sgn = sign::negative;
+  }
+  return res;
+}
+
+template <container C>
+constexpr z<C> root(const z<C>& num, int k) {
+  if (is_negative(num)) [[unlikely]] {
+    throw negative_radicand_error{};
+  }
+
+  if (k == 0) {
+    return details::one<C>();
+  } else if (k == 1 || is_zero(num)) {
+    return num;
+  }
+
+  z<C> x0 = num;  // guess a better intial value for faster convergence.
+  z<C> x1;
+  for (;;) {
+    auto t1 = mul_n(create<C>(k - 1), x0);
+    auto [t2, _] = div(num, pow(x0, k - 1));
+    x1 = div_n(add_n(t1, t2), create<C>(k)).q;
+    if (cmp_n(x1, x0) >= 0) {
+      return x1;
+    }
+    x0 = x1;
+  }
+}
+
 }  // namespace epx
 
 #endif

src/ut/chars_tests.cpp

@@ -18,25 +18,25 @@ TEST(chars_tests, try_from_decimal_chars) {
     auto opt_num = epx::try_from_chars<sz::container_type>("");
     EXPECT_TRUE(opt_num.has_value());
     EXPECT_TRUE(epx::is_zero(*opt_num));
-    EXPECT_TRUE(epx::is_positive(*opt_num));
+    EXPECT_FALSE(epx::is_negative(*opt_num));
   }
   {
     auto opt_num = epx::try_from_chars<sz::container_type>("0");
     EXPECT_TRUE(opt_num.has_value());
     EXPECT_TRUE(epx::is_zero(*opt_num));
-    EXPECT_TRUE(epx::is_positive(*opt_num));
+    EXPECT_FALSE(epx::is_negative(*opt_num));
   }
   {
     auto opt_num = epx::try_from_chars<sz::container_type>("-0");
     EXPECT_TRUE(opt_num.has_value());
     EXPECT_TRUE(epx::is_zero(*opt_num));
-    EXPECT_TRUE(epx::is_positive(*opt_num));
+    EXPECT_FALSE(epx::is_negative(*opt_num));
   }
   {
     auto opt_num = epx::try_from_chars<sz::container_type>("000000");
     EXPECT_TRUE(opt_num.has_value());
     EXPECT_TRUE(epx::is_zero(*opt_num));
-    EXPECT_TRUE(epx::is_positive(*opt_num));
+    EXPECT_FALSE(epx::is_negative(*opt_num));
   }
   {
     auto opt_num = epx::try_from_chars<sz::container_type>("1");

src/ut/def.hpp

@@ -4,6 +4,7 @@
 #pragma once
 
 // std
+#include <concepts>
 #include <cstdint>
 #include <string_view>
 #include <vector>
@@ -18,5 +19,9 @@ using mz = epx::z<std::vector<uint16_t>>;
 using lz = epx::z<std::vector<uint32_t>>;
 
 constexpr sz stosz(std::string_view chars) { return epx::try_from_chars<sz::container_type>(chars).value(); }
+template <std::integral T>
+constexpr sz create_sz(T val) {
+  return epx::create<sz::container_type>(val);
+}
 
 }  // namespace epxut

src/ut/r_tests.cpp

@@ -241,4 +241,36 @@ TEST(r_tests, inv) {
   }
 }
 
+TEST(r_tests, root) {
+  {
+    auto x = epx::make_q(stosz("2"), stosz("3"));
+    EXPECT_THROW(epx::root(x, 0).approx(2).get(), epx::kthroot_too_small_error);
+    EXPECT_THROW(epx::root(x, 1).approx(2).get(), epx::kthroot_too_small_error);
+  }
+  {
+    auto zero = epx::make_q(stosz("0"), stosz("1"));
+    auto expr = epx::root(zero, 2);
+    EXPECT_EQ("0.00000", epx::to_string(expr, 5));
+    EXPECT_EQ("0.0000000000", epx::to_string(expr, 10));
+  }
+  {
+    auto x = epx::make_q(stosz("1"), stosz("1"));
+    EXPECT_EQ("1.00000", epx::to_string(epx::root(x, 2), 5));
+    EXPECT_EQ("1.000", epx::to_string(epx::root(x, 11), 3));
+  }
+  {
+    auto x = epx::make_q(stosz("2"), stosz("1"));
+    auto expr = epx::root(x, 2);
+    EXPECT_EQ("1.41421", epx::to_string(expr, 5));
+    EXPECT_EQ("1.4142135624", epx::to_string(expr, 10));
+    EXPECT_EQ("1.4142135624", epx::to_string(expr, 10));
+    EXPECT_EQ("1.4142135623730950488016887242096980785697", epx::to_string(expr, 40));
+  }
+  {
+    auto x = epx::make_q(stosz("34012224"), stosz("1000000"));
+    auto expr = epx::root(x, 6);
+    EXPECT_EQ("1.8000000000", epx::to_string(expr, 10));
+  }
+}
+
 }  // namespace epxut

src/ut/z_tests.cpp

@@ -15,49 +15,66 @@
 
 namespace epxut {
 
+TEST(z_tests, create) {
+  {
+    auto num = epx::create<sz::container_type>(0);
+    EXPECT_TRUE(epx::is_zero(num));
+  }
+  {
+    auto num = epx::create<sz::container_type>(0x030201);
+    sz expected{.digits = {1, 2, 3}};
+    EXPECT_EQ(expected, num);
+  }
+  {
+    auto num = epx::create<mz::container_type>(-2);
+    mz expected{.digits = {2}, .sgn = epx::sign::negative};
+    EXPECT_EQ(expected, num);
+  }
+}
+
 TEST(z_tests, normalize) {
   {
     sz zero;
     EXPECT_TRUE(epx::is_zero(zero));
-    EXPECT_TRUE(epx::is_positive(zero));
+    EXPECT_FALSE(epx::is_negative(zero));
     epx::normalize(zero);
     EXPECT_TRUE(epx::is_zero(zero));
-    EXPECT_TRUE(epx::is_positive(zero));
+    EXPECT_FALSE(epx::is_negative(zero));
   }
   {
     mz zero;
     EXPECT_TRUE(epx::is_zero(zero));
-    EXPECT_TRUE(epx::is_positive(zero));
+    EXPECT_FALSE(epx::is_negative(zero));
 
     zero.digits = {0, 0};
     EXPECT_FALSE(epx::is_zero(zero));
-    EXPECT_TRUE(epx::is_positive(zero));
+    EXPECT_FALSE(epx::is_negative(zero));
     epx::normalize(zero);
     epx::add_n(zero, zero);
     EXPECT_TRUE(epx::is_zero(zero));
-    EXPECT_TRUE(epx::is_positive(zero));
+    EXPECT_FALSE(epx::is_negative(zero));
   }
   {
     lz zero;
     EXPECT_TRUE(epx::is_zero(zero));
-    EXPECT_TRUE(epx::is_positive(zero));
+    EXPECT_FALSE(epx::is_negative(zero));
 
     zero.digits = {0, 0};
     zero.sgn = epx::sign::negative;
     EXPECT_FALSE(epx::is_zero(zero));
-    EXPECT_FALSE(epx::is_positive(zero));
+    EXPECT_TRUE(epx::is_negative(zero));
     epx::normalize(zero);
     EXPECT_TRUE(epx::is_zero(zero));
-    EXPECT_TRUE(epx::is_positive(zero));
+    EXPECT_FALSE(epx::is_negative(zero));
   }
   {
     sz num = {.digits = {0, 1, 2, 0, 0}, .sgn = epx::sign::negative};
     EXPECT_FALSE(epx::is_zero(num));
-    EXPECT_FALSE(epx::is_positive(num));
+    EXPECT_TRUE(epx::is_negative(num));
     epx::normalize(num);
     auto expected_digits = sz::container_type{0, 1, 2};
     EXPECT_EQ(expected_digits, num.digits);
-    EXPECT_FALSE(epx::is_positive(num));
+    EXPECT_TRUE(epx::is_negative(num));
   }
 }
 
@@ -123,7 +140,7 @@ TEST(z_tests, mul) {
     EXPECT_TRUE(epx::is_zero(epx::mul(zero, zero)));
     EXPECT_TRUE(epx::is_zero(epx::mul(zero, one)));
     EXPECT_TRUE(epx::is_zero(epx::mul(one, zero)));
-    EXPECT_TRUE(epx::is_positive(epx::mul(minus_one, zero)));
+    EXPECT_FALSE(epx::is_negative(epx::mul(minus_one, zero)));
 
     EXPECT_EQ(one, epx::mul(one, one));
     EXPECT_EQ(one, epx::mul(minus_one, minus_one));
@@ -313,7 +330,7 @@ TEST(z_tests, mul_4exp) {
     epx::mul_4exp(num, -32);  // Large negative exponent
     sz expected{};
     EXPECT_EQ(expected, num);
-    EXPECT_TRUE(epx::is_positive(num));  // Should reset sign to positive
+    EXPECT_FALSE(epx::is_negative(num));  // Should reset sign to positive
   }
   {
     sz num{};
@@ -324,4 +341,60 @@ TEST(z_tests, mul_4exp) {
   }
 }
 
+TEST(z_tests, pow) {
+  {
+    sz zero;
+    sz one = {.digits = {1}};
+    auto pwr = epx::pow(zero, 0);  // 0^0 = 1
+    EXPECT_EQ(one, pwr);
+  }
+  {
+    sz zero;
+    auto pwr = epx::pow(zero, 1);  // 0^1 = 0
+    EXPECT_EQ(zero, pwr);
+  }
+  {
+    auto num = create_sz(2);
+    auto expected = create_sz(64);
+    auto pwr = epx::pow(num, 6);  //  2^6 = 64
+    EXPECT_EQ(expected, pwr);
+  }
+  {
+    sz one = {.digits = {1}};
+    sz minus_one{.digits = {1}, .sgn = epx::sign::negative};
+    EXPECT_EQ(minus_one, epx::pow(minus_one, 9));  // (-1)^9 = -1
+    EXPECT_EQ(one, epx::pow(minus_one, 10));       // (-1)^10 = 1
+  }
+}
+
+TEST(z_tests, root) {
+  {
+    sz zero;
+    EXPECT_TRUE(epx::is_zero(epx::root(zero, 1)));
+    EXPECT_TRUE(epx::is_zero(epx::root(zero, 2)));
+    EXPECT_TRUE(epx::is_zero(epx::root(zero, 7)));
+  }
+  {
+    sz one = {.digits = {1}};
+    EXPECT_EQ(one, epx::root(sz{}, 0));
+    EXPECT_EQ(one, epx::root(one, 0));
+    EXPECT_EQ(one, epx::root(create_sz(2), 0));
+    EXPECT_EQ(one, epx::root(create_sz(9), 0));
+  }
+  {
+    EXPECT_EQ(create_sz(2), epx::root(create_sz(4), 2));
+    EXPECT_EQ(create_sz(3), epx::root(create_sz(9), 2));
+    EXPECT_EQ(create_sz(19), epx::root(create_sz(361), 2));
+    EXPECT_EQ(create_sz(1000), epx::root(create_sz(1000000), 2));
+  }
+  {
+    EXPECT_EQ(create_sz(5), epx::root(create_sz(125), 3));
+    EXPECT_EQ(create_sz(37398), epx::root(create_sz(1956111062177043216), 4));
+  }
+  {
+    EXPECT_EQ(create_sz(111), epx::root(create_sz(12345), 2));
+    EXPECT_EQ(create_sz(23), epx::root(create_sz(12345), 3));
+  }
+}
+
 }  // namespace epxut