Add free function for calculation of simple continued fraction coeffs

Author: mborland

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 <limits>
 #include <stdexcept>
 #include <sstream>
+#include <utility>
+#include <cstdint>
 
 #include <boost/math/tools/is_standalone.hpp>
 #ifndef BOOST_MATH_STANDALONE
 #include <boost/config.hpp>
 #ifdef BOOST_NO_CXX17_IF_CONSTEXPR
-#error "The header <boost/math/norms.hpp> can only be used in C++17 and later."
+#error "The header <boost/math/simple_continued_fraction.hpp> 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<Real, 7> logs{std::numeric_limits<Real>::quiet_NaN(), Real(0), log(static_cast<Real>(2)), log(static_cast<Real>(3)), log(static_cast<Real>(4)), log(static_cast<Real>(5)), log(static_cast<Real>(6))};
+         const std::array<Real, 7> logs{std::numeric_limits<Real>::quiet_NaN(), static_cast<Real>(0), log(static_cast<Real>(2)), log(static_cast<Real>(3)), log(static_cast<Real>(4)), log(static_cast<Real>(5)), log(static_cast<Real>(6))};
          Real log_prod = 0;
          for (size_t i = 1; i < b_.size(); ++i) {
             if (b_[i] < static_cast<Z>(logs.size())) {
@@ -138,6 +141,11 @@ class simple_continued_fraction {
     const std::vector<Z>& partial_denominators() const {
       return b_;
     }
+
+    inline std::vector<Z>&& get_data() noexcept
+    {
+        return std::move(b_);
+    }
 
     template<typename T, typename Z2>
     friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction<T, Z2>& scf);
@@ -171,6 +179,12 @@ std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>&
    return out;
 }
 
+template<typename Real, typename Z = std::int64_t>
+inline auto simple_continued_fraction_coefficients(Real x)
+{
+    auto temp = simple_continued_fraction(x);
+    return temp.get_data();
+}
 
 }
 #endif

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 <typename Real>
+void test_git_issue_970()
+{
+    using boost::math::tools::simple_continued_fraction_coefficients;
+
+    auto coefs1 = simple_continued_fraction_coefficients(static_cast<Real>(3.14155L)); // [3; 7, 15, 2, 7, 1, 4, 2]
+    auto coefs2 = simple_continued_fraction_coefficients(static_cast<Real>(3.14165L)); // [3; 7, 16, 1, 3, 4, 2, 4]
+
+    const std::size_t max_size = (std::min)(coefs1.size(), coefs2.size());
+    std::vector<std::int64_t> 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<float128>();
     test_khinchin<float128>();
     #endif
+
+    test_git_issue_970<float>();
+    test_git_issue_970<double>();
+    test_git_issue_970<long double>();
+
     return boost::math::test::report_errors();
 }