From 08f04908752a3263566a29513256190425a43ec3 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 4 Apr 2023 14:58:15 +0200 Subject: [PATCH 1/6] Add free function for calculation of simple continued fraction coeffs --- .../math/tools/simple_continued_fraction.hpp | 18 ++++++++- test/simple_continued_fraction_test.cpp | 38 +++++++++++++++++++ 2 files changed, 54 insertions(+), 2 deletions(-) diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp index 0fe17fc59d..f904c5f279 100644 --- a/include/boost/math/tools/simple_continued_fraction.hpp +++ b/include/boost/math/tools/simple_continued_fraction.hpp @@ -1,4 +1,5 @@ // (C) Copyright Nick Thompson 2020. +// (C) Copyright Matt Borland 2023. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -14,12 +15,14 @@ #include #include #include +#include +#include #include #ifndef BOOST_MATH_STANDALONE #include #ifdef BOOST_NO_CXX17_IF_CONSTEXPR -#error "The header can only be used in C++17 and later." +#error "The header can only be used in C++17 and later." #endif #endif @@ -108,7 +111,7 @@ class simple_continued_fraction { // Precompute the most probable logarithms. See the Gauss-Kuzmin distribution for details. // Example: b_i = 1 has probability -log_2(3/4) ~ .415: // A random partial denominator has ~80% chance of being in this table: - const std::array logs{std::numeric_limits::quiet_NaN(), Real(0), log(static_cast(2)), log(static_cast(3)), log(static_cast(4)), log(static_cast(5)), log(static_cast(6))}; + 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) { if (b_[i] < static_cast(logs.size())) { @@ -138,6 +141,11 @@ class simple_continued_fraction { const std::vector& partial_denominators() const { return b_; } + + inline std::vector&& get_data() noexcept + { + return std::move(b_); + } template friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf); @@ -171,6 +179,12 @@ std::ostream& operator<<(std::ostream& out, simple_continued_fraction& return out; } +template +inline auto simple_continued_fraction_coefficients(Real x) +{ + auto temp = simple_continued_fraction(x); + return temp.get_data(); +} } #endif diff --git a/test/simple_continued_fraction_test.cpp b/test/simple_continued_fraction_test.cpp index fedb4e5a30..9f06aad2cf 100644 --- a/test/simple_continued_fraction_test.cpp +++ b/test/simple_continued_fraction_test.cpp @@ -1,5 +1,6 @@ /* * Copyright Nick Thompson, 2020 + * Copyright Matt Borland, 2023 * Use, modification and distribution are subject to the * Boost Software License, Version 1.0. (See accompanying file * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -131,6 +132,38 @@ void test_khinchin() CHECK_ULP_CLOSE(Real(2), Km1, 10); } +template +void test_git_issue_970() +{ + using boost::math::tools::simple_continued_fraction_coefficients; + + auto coefs1 = simple_continued_fraction_coefficients(static_cast(3.14155L)); // [3; 7, 15, 2, 7, 1, 4, 2] + auto coefs2 = simple_continued_fraction_coefficients(static_cast(3.14165L)); // [3; 7, 16, 1, 3, 4, 2, 4] + + const std::size_t max_size = (std::min)(coefs1.size(), coefs2.size()); + std::vector coefs; + coefs.reserve(max_size); + + for (std::size_t i = 0; i < max_size; ++i) + { + const auto c1 = coefs1[i]; + const auto c2 = coefs2[i]; + if (c1 == c2) + { + coefs.emplace_back(c1); + continue; + } + + coefs.emplace_back((std::min)(c1, c2) + 1); + break; + } + + // Result is [3; 7, 16] + CHECK_EQUAL(coefs.size(), 3UL); + CHECK_EQUAL(coefs[0], INT64_C(3)); + CHECK_EQUAL(coefs[1], INT64_C(7)); + CHECK_EQUAL(coefs[2], INT64_C(16)); +} int main() { @@ -157,5 +190,10 @@ int main() test_simple(); test_khinchin(); #endif + + test_git_issue_970(); + test_git_issue_970(); + test_git_issue_970(); + return boost::math::test::report_errors(); } From b547ac7ed05bbeb134911685003a232f76ca0c62 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 4 Apr 2023 15:09:27 +0200 Subject: [PATCH 2/6] Update docs --- doc/internals/simple_continued_fraction.qbk | 42 ++++++++++++++++++++- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/doc/internals/simple_continued_fraction.qbk b/doc/internals/simple_continued_fraction.qbk index 45896ffd74..e99b6eaefe 100644 --- a/doc/internals/simple_continued_fraction.qbk +++ b/doc/internals/simple_continued_fraction.qbk @@ -1,5 +1,6 @@ [/ Copyright Nick Thompson, 2020 + Copyright Matt Borland, 2023 Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt). @@ -19,9 +20,17 @@ Real khinchin_harmonic_mean() const; - template - friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf); + const std::vector& partial_denominators() const; + + inline std::vector&& get_data() noexcept; + + template + friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction& scf); }; + + template + inline std::vector simple_continued_fraction_coefficients(Real x); + } @@ -47,6 +56,35 @@ This is because when examining known values like π, it creates a large number o It may be the case the a few incorrect partial convergents is harmless, but we compute continued fractions because we would like to do something with them. One sensible thing to do it to ask whether the number is in some sense "random"; a question that can be partially answered by computing the Khinchin geometric mean +If you only require the coefficients of the simple continued fraction for example in the calculation of [@https://en.wikipedia.org/wiki/Continued_fraction#Best_rational_approximations best rational approximations] there is a free function for that. + +An example of this calculation follows: + + using boost::math::tools::simple_continued_fraction_coefficients; + + auto coefs1 = simple_continued_fraction_coefficients(static_cast(3.14155L)); // [3; 7, 15, 2, 7, 1, 4, 2] + auto coefs2 = simple_continued_fraction_coefficients(static_cast(3.14165L)); // [3; 7, 16, 1, 3, 4, 2, 4] + + const std::size_t max_size = (std::min)(coefs1.size(), coefs2.size()); + std::vector coefs; + coefs.reserve(max_size); + + for (std::size_t i = 0; i < max_size; ++i) + { + const auto c1 = coefs1[i]; + const auto c2 = coefs2[i]; + if (c1 == c2) + { + coefs.emplace_back(c1); + continue; + } + + coefs.emplace_back((std::min)(c1, c2) + 1); + break; + } + + // Result is [3; 7, 16] + [$../equations/khinchin_geometric.svg] and Khinchin harmonic mean From 1f3a079945476fbf57e9548c03e58ad120ccecc0 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 4 Apr 2023 16:16:57 +0200 Subject: [PATCH 3/6] Fix casting for MSVC --- test/simple_continued_fraction_test.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/simple_continued_fraction_test.cpp b/test/simple_continued_fraction_test.cpp index 9f06aad2cf..65afef59c7 100644 --- a/test/simple_continued_fraction_test.cpp +++ b/test/simple_continued_fraction_test.cpp @@ -159,10 +159,9 @@ void test_git_issue_970() } // Result is [3; 7, 16] - CHECK_EQUAL(coefs.size(), 3UL); - CHECK_EQUAL(coefs[0], INT64_C(3)); - CHECK_EQUAL(coefs[1], INT64_C(7)); - CHECK_EQUAL(coefs[2], INT64_C(16)); + CHECK_EQUAL(coefs[0], static_cast(3)); + CHECK_EQUAL(coefs[1], static_cast(7)); + CHECK_EQUAL(coefs[2], static_cast(16)); } int main() 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 4/6] 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 5/6] 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 6/6] 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