From d83bd538381eb8f047dd2ab49be6c7bc4beaacb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Kol=C3=A1rik?= Date: Tue, 11 Apr 2023 16:12:28 +0200 Subject: [PATCH 1/3] Typedef for Z `int_type`; replaced member variable `x_` with `precision_` --- .../math/tools/simple_continued_fraction.hpp | 71 ++++++++++--------- 1 file changed, 39 insertions(+), 32 deletions(-) diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp index f904c5f279..e32b1bfaf2 100644 --- a/include/boost/math/tools/simple_continued_fraction.hpp +++ b/include/boost/math/tools/simple_continued_fraction.hpp @@ -35,14 +35,25 @@ namespace boost::math::tools { template class simple_continued_fraction { public: - simple_continued_fraction(Real x) : x_{x} { + typedef Z int_type; + + simple_continued_fraction(Real x) { using std::floor; using std::abs; using std::sqrt; using std::isfinite; + const Real orig_x = x; if (!isfinite(x)) { - throw std::domain_error("Cannot convert non-finites into continued fractions."); + throw std::domain_error("Cannot convert non-finites into continued fractions."); + } + + constexpr int p = std::numeric_limits::max_digits10; + if constexpr (p == 2147483647) { + precision_ = orig_x.backend().precision(); + } else { + precision_ = p; } + b_.reserve(50); Real bj = floor(x); b_.push_back(static_cast(bj)); @@ -57,12 +68,11 @@ class simple_continued_fraction { } Real C = f; Real D = 0; - int i = 0; - // the "1 + i++" lets the error bound grow slowly with the number of convergents. + // the "1 + i" lets the error bound grow slowly with the number of convergents. // I have not worked out the error propagation of the Modified Lentz's method to see if it does indeed grow at this rate. // Numerical Recipes claims that no one has worked out the error analysis of the modified Lentz's method. - while (abs(f - x_) >= (1 + i++)*std::numeric_limits::epsilon()*abs(x_)) - { + const Real eps_abs_orig_x = std::numeric_limits::epsilon()*abs(orig_x); + for (int i = 0; abs(f - orig_x) >= (1 + i)*eps_abs_orig_x; ++i) { bj = floor(x); b_.push_back(static_cast(bj)); x = 1/(x-bj); @@ -81,12 +91,13 @@ class simple_continued_fraction { // The shorter representation is considered the canonical representation, // so if we compute a non-canonical representation, change it to canonical: if (b_.size() > 2 && b_.back() == 1) { - b_[b_.size() - 2] += 1; - b_.resize(b_.size() - 1); + b_.pop_back(); + b_.back() += 1; } b_.shrink_to_fit(); - - for (size_t i = 1; i < b_.size(); ++i) { + + const size_t size_ = b_.size(); + for (size_t i = 1; i < size_; ++i) { if (b_[i] <= 0) { std::ostringstream oss; oss << "Found a negative partial denominator: b[" << i << "] = " << b_[i] << "." @@ -101,9 +112,10 @@ class simple_continued_fraction { } } } - + Real khinchin_geometric_mean() const { - if (b_.size() == 1) { + const size_t size_ = b_.size(); + if (size_ == 1) { return std::numeric_limits::quiet_NaN(); } using std::log; @@ -113,7 +125,7 @@ class simple_continued_fraction { // A random partial denominator has ~80% chance of being in this table: const std::array logs{std::numeric_limits::quiet_NaN(), static_cast(0), log(static_cast(2)), log(static_cast(3)), log(static_cast(4)), log(static_cast(5)), log(static_cast(6))}; Real log_prod = 0; - for (size_t i = 1; i < b_.size(); ++i) { + for (size_t i = 1; i < size_; ++i) { if (b_[i] < static_cast(logs.size())) { log_prod += logs[b_[i]]; } @@ -122,49 +134,44 @@ class simple_continued_fraction { log_prod += log(static_cast(b_[i])); } } - log_prod /= (b_.size()-1); + log_prod /= (size_-1); return exp(log_prod); } - + Real khinchin_harmonic_mean() const { - if (b_.size() == 1) { + const size_t size_ = b_.size(); + if (size_ == 1) { return std::numeric_limits::quiet_NaN(); } - Real n = b_.size() - 1; + Real n = size_ - 1; Real denom = 0; - for (size_t i = 1; i < b_.size(); ++i) { + for (size_t i = 1; i < size_; ++i) { denom += 1/static_cast(b_[i]); } return n/denom; } - + + // Note that this also includes the integer part (i.e. all the coefficients) const std::vector& partial_denominators() const { return b_; } - inline std::vector&& get_data() noexcept - { + inline std::vector&& get_data() noexcept { return std::move(b_); } - + template friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf); - private: - const Real x_; - std::vector b_; + std::vector b_{}; + + int precision_{}; }; template std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf) { - constexpr const int p = std::numeric_limits::max_digits10; - if constexpr (p == 2147483647) { - out << std::setprecision(scf.x_.backend().precision()); - } else { - out << std::setprecision(p); - } - + out << std::setprecision(scf.precision_); out << "[" << scf.b_.front(); if (scf.b_.size() > 1) { From 25ab7c57c889fe3194198fc1d7a0c145ff2932a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Kol=C3=A1rik?= Date: Tue, 11 Apr 2023 17:03:59 +0200 Subject: [PATCH 2/3] Implemented `simple_continued_fraction(std::vector)` --- .../math/tools/simple_continued_fraction.hpp | 57 +++++++++++++------ 1 file changed, 39 insertions(+), 18 deletions(-) diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp index e32b1bfaf2..2f7541d25b 100644 --- a/include/boost/math/tools/simple_continued_fraction.hpp +++ b/include/boost/math/tools/simple_continued_fraction.hpp @@ -37,21 +37,34 @@ class simple_continued_fraction { public: typedef Z int_type; - simple_continued_fraction(Real x) { + simple_continued_fraction(std::vector data) : b_{std::move(data)} { + const size_t size_ = b_.size(); + if (size_ == 0) { + throw std::length_error("Array of coefficients is empty."); + } + + for (size_t i = 1; i < size_; ++i) { + if (b_[i] <= 0) { + std::ostringstream oss; + oss << "Found a negative partial denominator: b[" << i << "] = " << b_[i] << "."; + throw std::domain_error(oss.str()); + } + } + + canonicalize(); + } + + simple_continued_fraction(Real x) : b_{} { using std::floor; using std::abs; using std::sqrt; using std::isfinite; - const Real orig_x = x; if (!isfinite(x)) { throw std::domain_error("Cannot convert non-finites into continued fractions."); } - constexpr int p = std::numeric_limits::max_digits10; - if constexpr (p == 2147483647) { - precision_ = orig_x.backend().precision(); - } else { - precision_ = p; + if constexpr (std_precision == 2147483647) { + precision_ = x.backend().precision(); } b_.reserve(50); @@ -61,6 +74,8 @@ class simple_continued_fraction { b_.shrink_to_fit(); return; } + + const Real orig_x = x; x = 1/(x-bj); Real f = bj; if (bj == 0) { @@ -87,14 +102,7 @@ class simple_continued_fraction { D = 1/D; f *= (C*D); } - // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1]. - // The shorter representation is considered the canonical representation, - // so if we compute a non-canonical representation, change it to canonical: - if (b_.size() > 2 && b_.back() == 1) { - b_.pop_back(); - b_.back() += 1; - } - b_.shrink_to_fit(); + canonicalize(); const size_t size_ = b_.size(); for (size_t i = 1; i < size_; ++i) { @@ -163,9 +171,22 @@ class simple_continued_fraction { template friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf); private: - std::vector b_{}; + static constexpr int std_precision = std::numeric_limits::max_digits10; + + void canonicalize() { + // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1]. + // The shorter representation is considered the canonical representation, + // so if we compute a non-canonical representation, change it to canonical: + if (b_.size() > 2 && b_.back() == 1) { + b_.pop_back(); + b_.back() += 1; + } + b_.shrink_to_fit(); + } + + std::vector b_; - int precision_{}; + int precision_{std_precision}; }; @@ -189,7 +210,7 @@ std::ostream& operator<<(std::ostream& out, simple_continued_fraction& template inline auto simple_continued_fraction_coefficients(Real x) { - auto temp = simple_continued_fraction(x); + auto temp = simple_continued_fraction(x); return temp.get_data(); } From c4c4592d07d297ba804451a9b6082518e3153260 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Kol=C3=A1rik?= Date: Tue, 11 Apr 2023 17:30:19 +0200 Subject: [PATCH 3/3] Implemented templated free function `Rational to_rational(simple_continued_fraction)` --- .../math/tools/simple_continued_fraction.hpp | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp index 2f7541d25b..5402d2f6ff 100644 --- a/include/boost/math/tools/simple_continued_fraction.hpp +++ b/include/boost/math/tools/simple_continued_fraction.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #ifndef BOOST_MATH_STANDALONE @@ -214,5 +215,34 @@ inline auto simple_continued_fraction_coefficients(Real x) return temp.get_data(); } +// Can be used with `boost::rational` from +template +inline Rational to_rational(const simple_continued_fraction& scf) +{ + using int_t = typename Rational::int_type; + + auto& coefs = scf.partial_denominators(); + const size_t size_ = coefs.size(); + assert(size_ >= 1); + if (size_ == 1) return static_cast(coefs[0]); + + // p0 = a0, p1 = a1.a0 + 1, pn = an.pn-1 + pn-2 for 2 <= n + // q0 = 1, q1 = a1, qn = an.qn-1 + qn-2 for 2 <= n + + int_t p0 = coefs[0]; + int_t p1 = p0*coefs[1] + 1; + int_t q0 = 1; + int_t q1 = coefs[1]; + for (size_t i = 2; i < size_; ++i) { + const Z cn = coefs[i]; + const int_t pn = cn*p1 + p0; + const int_t qn = cn*q1 + q0; + p0 = std::exchange(p1, pn); + q0 = std::exchange(q1, qn); + } + + return {p1, q1}; +} + } #endif