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 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..65afef59c7 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,37 @@ 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[0], static_cast(3)); + CHECK_EQUAL(coefs[1], static_cast(7)); + CHECK_EQUAL(coefs[2], static_cast(16)); +} int main() { @@ -157,5 +189,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(); }