Skip to content

Commit f1dabc5

Browse files
tian-lt-personaltian-lt-personal
authored
algo.6 kth root. (#19)
1 parent 22b44d6 commit f1dabc5

File tree

8 files changed

+224
-21
lines changed

8 files changed

+224
-21
lines changed

‎src/epsilon/chars.hpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -107,9 +107,9 @@ constexpr std::string to_string(r<C> num, unsigned int k) {
107107
d.sgn = sgn; // restore the sign
108108

109109
auto s = to_string(d);
110-
size_t len = is_positive(d) ? s.length() : s.length() - 1;
110+
size_t len = is_negative(d) ? s.length() - 1 : s.length();
111111
if (len <= k) {
112-
s.insert(is_positive(d) ? 0 : 1, k - len + 1, '0');
112+
s.insert(is_negative(d) ? 1 : 0, k - len + 1, '0');
113113
}
114114
if (k > 0) {
115115
s.insert(s.length() - k, 1, '.');

‎src/epsilon/r.hpp‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,21 @@ constexpr r<C> inv(r<C> x) {
148148
}};
149149
}
150150

151+
template <container C>
152+
constexpr r<C> root(r<C> x, int k) {
153+
return r<C>{[x = std::move(x), k](int n) -> coro::lazy<z<C>> {
154+
if (k < 2) [[unlikely]] {
155+
throw kthroot_too_small_error{};
156+
}
157+
158+
auto xn = co_await x.approx(n * k);
159+
if (is_negative(xn)) [[unlikely]] {
160+
throw negative_radicand_error{};
161+
}
162+
co_return root(xn, k);
163+
}};
164+
}
165+
151166
} // namespace epx
152167

153168
#endif // EPSILON_INC_R_HPP

‎src/epsilon/t.hpp‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,18 @@ struct precision_overflow_error : public std::runtime_error {
6161
precision_overflow_error() : std::runtime_error("epx: precision overflow error") {}
6262
};
6363

64+
struct negative_radicand_error : public std::runtime_error {
65+
negative_radicand_error() : std::runtime_error("epx: negative radicand error") {}
66+
};
67+
68+
struct negative_zpow_error : public std::runtime_error {
69+
negative_zpow_error() : std::runtime_error("epx: negative power for z error") {}
70+
};
71+
72+
struct kthroot_too_small_error : public std::runtime_error {
73+
kthroot_too_small_error() : std::runtime_error("epx: kth-root too small error") {}
74+
};
75+
6476
} // namespace epx
6577

6678
#endif

‎src/epsilon/z.hpp‎

Lines changed: 69 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -117,14 +117,36 @@ constexpr auto bit_shift(C& digits, int offset) {
117117
}
118118
} // namespace details
119119

120+
template <container C, std::integral T>
121+
constexpr z<C> create(T val) {
122+
z<C> res;
123+
if (val == 0) {
124+
return res;
125+
}
126+
if constexpr (std::is_signed_v<T>) {
127+
res.sgn = val < 0 ? sign::negative : sign::positive;
128+
val = val < 0 ? -val : val;
129+
}
130+
using D = typename z<C>::digit_type;
131+
if constexpr (sizeof(T) <= sizeof(D)) {
132+
res.digits.push_back(static_cast<D>(val));
133+
} else {
134+
while (val > 0) {
135+
res.digits.push_back(static_cast<D>(val));
136+
val >>= (sizeof(D) * CHAR_BIT);
137+
}
138+
}
139+
return res;
140+
}
141+
120142
template <container C>
121143
constexpr bool is_zero(const z<C>& num) noexcept {
122144
return std::ranges::size(num.digits) == 0;
123145
}
124146

125147
template <container C>
126-
constexpr bool is_positive(const z<C>& num) noexcept {
127-
return num.sgn == sign::positive;
148+
constexpr bool is_negative(const z<C>& num) noexcept {
149+
return num.sgn == sign::negative;
128150
}
129151

130152
template <container C>
@@ -141,7 +163,7 @@ constexpr z<C>& negate(z<C>& num) noexcept {
141163
if (is_zero(num)) {
142164
return num;
143165
}
144-
num.sgn = is_positive(num) ? sign::negative : sign::positive;
166+
num.sgn = is_negative(num) ? sign::positive : sign::negative;
145167
return num;
146168
}
147169

@@ -503,6 +525,50 @@ constexpr z<C> mul_4exp(const z<C>& val, int exp) {
503525
return mul_4exp(num, exp);
504526
}
505527

528+
template <container C>
529+
constexpr z<C> pow(const z<C>& num, int exp) {
530+
assert(exp >= 0);
531+
z<C> res = details::one<C>();
532+
if (exp > 0) {
533+
for (int i = 0; i < exp; ++i) {
534+
res = mul_n(res, num);
535+
}
536+
} else if (exp == 0) {
537+
return res;
538+
} else [[unlikely]] {
539+
throw negative_zpow_error{};
540+
}
541+
if (is_negative(num) && (exp % 2 == 1)) {
542+
res.sgn = sign::negative;
543+
}
544+
return res;
545+
}
546+
547+
template <container C>
548+
constexpr z<C> root(const z<C>& num, int k) {
549+
if (is_negative(num)) [[unlikely]] {
550+
throw negative_radicand_error{};
551+
}
552+
553+
if (k == 0) {
554+
return details::one<C>();
555+
} else if (k == 1 || is_zero(num)) {
556+
return num;
557+
}
558+
559+
z<C> x0 = num; // guess a better intial value for faster convergence.
560+
z<C> x1;
561+
for (;;) {
562+
auto t1 = mul_n(create<C>(k - 1), x0);
563+
auto [t2, _] = div(num, pow(x0, k - 1));
564+
x1 = div_n(add_n(t1, t2), create<C>(k)).q;
565+
if (cmp_n(x1, x0) >= 0) {
566+
return x1;
567+
}
568+
x0 = x1;
569+
}
570+
}
571+
506572
} // namespace epx
507573

508574
#endif

‎src/ut/chars_tests.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -18,25 +18,25 @@ TEST(chars_tests, try_from_decimal_chars) {
1818
auto opt_num = epx::try_from_chars<sz::container_type>("");
1919
EXPECT_TRUE(opt_num.has_value());
2020
EXPECT_TRUE(epx::is_zero(*opt_num));
21-
EXPECT_TRUE(epx::is_positive(*opt_num));
21+
EXPECT_FALSE(epx::is_negative(*opt_num));
2222
}
2323
{
2424
auto opt_num = epx::try_from_chars<sz::container_type>("0");
2525
EXPECT_TRUE(opt_num.has_value());
2626
EXPECT_TRUE(epx::is_zero(*opt_num));
27-
EXPECT_TRUE(epx::is_positive(*opt_num));
27+
EXPECT_FALSE(epx::is_negative(*opt_num));
2828
}
2929
{
3030
auto opt_num = epx::try_from_chars<sz::container_type>("-0");
3131
EXPECT_TRUE(opt_num.has_value());
3232
EXPECT_TRUE(epx::is_zero(*opt_num));
33-
EXPECT_TRUE(epx::is_positive(*opt_num));
33+
EXPECT_FALSE(epx::is_negative(*opt_num));
3434
}
3535
{
3636
auto opt_num = epx::try_from_chars<sz::container_type>("000000");
3737
EXPECT_TRUE(opt_num.has_value());
3838
EXPECT_TRUE(epx::is_zero(*opt_num));
39-
EXPECT_TRUE(epx::is_positive(*opt_num));
39+
EXPECT_FALSE(epx::is_negative(*opt_num));
4040
}
4141
{
4242
auto opt_num = epx::try_from_chars<sz::container_type>("1");

‎src/ut/def.hpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#pragma once
55

66
// std
7+
#include <concepts>
78
#include <cstdint>
89
#include <string_view>
910
#include <vector>
@@ -18,5 +19,9 @@ using mz = epx::z<std::vector<uint16_t>>;
1819
using lz = epx::z<std::vector<uint32_t>>;
1920

2021
constexpr sz stosz(std::string_view chars) { return epx::try_from_chars<sz::container_type>(chars).value(); }
22+
template <std::integral T>
23+
constexpr sz create_sz(T val) {
24+
return epx::create<sz::container_type>(val);
25+
}
2126

2227
} // namespace epxut

‎src/ut/r_tests.cpp‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -241,4 +241,36 @@ TEST(r_tests, inv) {
241241
}
242242
}
243243

244+
TEST(r_tests, root) {
245+
{
246+
auto x = epx::make_q(stosz("2"), stosz("3"));
247+
EXPECT_THROW(epx::root(x, 0).approx(2).get(), epx::kthroot_too_small_error);
248+
EXPECT_THROW(epx::root(x, 1).approx(2).get(), epx::kthroot_too_small_error);
249+
}
250+
{
251+
auto zero = epx::make_q(stosz("0"), stosz("1"));
252+
auto expr = epx::root(zero, 2);
253+
EXPECT_EQ("0.00000", epx::to_string(expr, 5));
254+
EXPECT_EQ("0.0000000000", epx::to_string(expr, 10));
255+
}
256+
{
257+
auto x = epx::make_q(stosz("1"), stosz("1"));
258+
EXPECT_EQ("1.00000", epx::to_string(epx::root(x, 2), 5));
259+
EXPECT_EQ("1.000", epx::to_string(epx::root(x, 11), 3));
260+
}
261+
{
262+
auto x = epx::make_q(stosz("2"), stosz("1"));
263+
auto expr = epx::root(x, 2);
264+
EXPECT_EQ("1.41421", epx::to_string(expr, 5));
265+
EXPECT_EQ("1.4142135624", epx::to_string(expr, 10));
266+
EXPECT_EQ("1.4142135624", epx::to_string(expr, 10));
267+
EXPECT_EQ("1.4142135623730950488016887242096980785697", epx::to_string(expr, 40));
268+
}
269+
{
270+
auto x = epx::make_q(stosz("34012224"), stosz("1000000"));
271+
auto expr = epx::root(x, 6);
272+
EXPECT_EQ("1.8000000000", epx::to_string(expr, 10));
273+
}
274+
}
275+
244276
} // namespace epxut

‎src/ut/z_tests.cpp‎

Lines changed: 85 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -15,49 +15,66 @@
1515

1616
namespace epxut {
1717

18+
TEST(z_tests, create) {
19+
{
20+
auto num = epx::create<sz::container_type>(0);
21+
EXPECT_TRUE(epx::is_zero(num));
22+
}
23+
{
24+
auto num = epx::create<sz::container_type>(0x030201);
25+
sz expected{.digits = {1, 2, 3}};
26+
EXPECT_EQ(expected, num);
27+
}
28+
{
29+
auto num = epx::create<mz::container_type>(-2);
30+
mz expected{.digits = {2}, .sgn = epx::sign::negative};
31+
EXPECT_EQ(expected, num);
32+
}
33+
}
34+
1835
TEST(z_tests, normalize) {
1936
{
2037
sz zero;
2138
EXPECT_TRUE(epx::is_zero(zero));
22-
EXPECT_TRUE(epx::is_positive(zero));
39+
EXPECT_FALSE(epx::is_negative(zero));
2340
epx::normalize(zero);
2441
EXPECT_TRUE(epx::is_zero(zero));
25-
EXPECT_TRUE(epx::is_positive(zero));
42+
EXPECT_FALSE(epx::is_negative(zero));
2643
}
2744
{
2845
mz zero;
2946
EXPECT_TRUE(epx::is_zero(zero));
30-
EXPECT_TRUE(epx::is_positive(zero));
47+
EXPECT_FALSE(epx::is_negative(zero));
3148

3249
zero.digits = {0, 0};
3350
EXPECT_FALSE(epx::is_zero(zero));
34-
EXPECT_TRUE(epx::is_positive(zero));
51+
EXPECT_FALSE(epx::is_negative(zero));
3552
epx::normalize(zero);
3653
epx::add_n(zero, zero);
3754
EXPECT_TRUE(epx::is_zero(zero));
38-
EXPECT_TRUE(epx::is_positive(zero));
55+
EXPECT_FALSE(epx::is_negative(zero));
3956
}
4057
{
4158
lz zero;
4259
EXPECT_TRUE(epx::is_zero(zero));
43-
EXPECT_TRUE(epx::is_positive(zero));
60+
EXPECT_FALSE(epx::is_negative(zero));
4461

4562
zero.digits = {0, 0};
4663
zero.sgn = epx::sign::negative;
4764
EXPECT_FALSE(epx::is_zero(zero));
48-
EXPECT_FALSE(epx::is_positive(zero));
65+
EXPECT_TRUE(epx::is_negative(zero));
4966
epx::normalize(zero);
5067
EXPECT_TRUE(epx::is_zero(zero));
51-
EXPECT_TRUE(epx::is_positive(zero));
68+
EXPECT_FALSE(epx::is_negative(zero));
5269
}
5370
{
5471
sz num = {.digits = {0, 1, 2, 0, 0}, .sgn = epx::sign::negative};
5572
EXPECT_FALSE(epx::is_zero(num));
56-
EXPECT_FALSE(epx::is_positive(num));
73+
EXPECT_TRUE(epx::is_negative(num));
5774
epx::normalize(num);
5875
auto expected_digits = sz::container_type{0, 1, 2};
5976
EXPECT_EQ(expected_digits, num.digits);
60-
EXPECT_FALSE(epx::is_positive(num));
77+
EXPECT_TRUE(epx::is_negative(num));
6178
}
6279
}
6380

@@ -123,7 +140,7 @@ TEST(z_tests, mul) {
123140
EXPECT_TRUE(epx::is_zero(epx::mul(zero, zero)));
124141
EXPECT_TRUE(epx::is_zero(epx::mul(zero, one)));
125142
EXPECT_TRUE(epx::is_zero(epx::mul(one, zero)));
126-
EXPECT_TRUE(epx::is_positive(epx::mul(minus_one, zero)));
143+
EXPECT_FALSE(epx::is_negative(epx::mul(minus_one, zero)));
127144

128145
EXPECT_EQ(one, epx::mul(one, one));
129146
EXPECT_EQ(one, epx::mul(minus_one, minus_one));
@@ -313,7 +330,7 @@ TEST(z_tests, mul_4exp) {
313330
epx::mul_4exp(num, -32); // Large negative exponent
314331
sz expected{};
315332
EXPECT_EQ(expected, num);
316-
EXPECT_TRUE(epx::is_positive(num)); // Should reset sign to positive
333+
EXPECT_FALSE(epx::is_negative(num)); // Should reset sign to positive
317334
}
318335
{
319336
sz num{};
@@ -324,4 +341,60 @@ TEST(z_tests, mul_4exp) {
324341
}
325342
}
326343

344+
TEST(z_tests, pow) {
345+
{
346+
sz zero;
347+
sz one = {.digits = {1}};
348+
auto pwr = epx::pow(zero, 0); // 0^0 = 1
349+
EXPECT_EQ(one, pwr);
350+
}
351+
{
352+
sz zero;
353+
auto pwr = epx::pow(zero, 1); // 0^1 = 0
354+
EXPECT_EQ(zero, pwr);
355+
}
356+
{
357+
auto num = create_sz(2);
358+
auto expected = create_sz(64);
359+
auto pwr = epx::pow(num, 6); // 2^6 = 64
360+
EXPECT_EQ(expected, pwr);
361+
}
362+
{
363+
sz one = {.digits = {1}};
364+
sz minus_one{.digits = {1}, .sgn = epx::sign::negative};
365+
EXPECT_EQ(minus_one, epx::pow(minus_one, 9)); // (-1)^9 = -1
366+
EXPECT_EQ(one, epx::pow(minus_one, 10)); // (-1)^10 = 1
367+
}
368+
}
369+
370+
TEST(z_tests, root) {
371+
{
372+
sz zero;
373+
EXPECT_TRUE(epx::is_zero(epx::root(zero, 1)));
374+
EXPECT_TRUE(epx::is_zero(epx::root(zero, 2)));
375+
EXPECT_TRUE(epx::is_zero(epx::root(zero, 7)));
376+
}
377+
{
378+
sz one = {.digits = {1}};
379+
EXPECT_EQ(one, epx::root(sz{}, 0));
380+
EXPECT_EQ(one, epx::root(one, 0));
381+
EXPECT_EQ(one, epx::root(create_sz(2), 0));
382+
EXPECT_EQ(one, epx::root(create_sz(9), 0));
383+
}
384+
{
385+
EXPECT_EQ(create_sz(2), epx::root(create_sz(4), 2));
386+
EXPECT_EQ(create_sz(3), epx::root(create_sz(9), 2));
387+
EXPECT_EQ(create_sz(19), epx::root(create_sz(361), 2));
388+
EXPECT_EQ(create_sz(1000), epx::root(create_sz(1000000), 2));
389+
}
390+
{
391+
EXPECT_EQ(create_sz(5), epx::root(create_sz(125), 3));
392+
EXPECT_EQ(create_sz(37398), epx::root(create_sz(1956111062177043216), 4));
393+
}
394+
{
395+
EXPECT_EQ(create_sz(111), epx::root(create_sz(12345), 2));
396+
EXPECT_EQ(create_sz(23), epx::root(create_sz(12345), 3));
397+
}
398+
}
399+
327400
} // namespace epxut

0 commit comments

Comments
 (0)