diff --git a/include/boost/math/tools/simple_continued_fraction.hpp b/include/boost/math/tools/simple_continued_fraction.hpp
index f904c5f279..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 <sstream>
 #include <utility>
 #include <cstdint>
+#include <cassert>
 
 #include <boost/math/tools/is_standalone.hpp>
 #ifndef BOOST_MATH_STANDALONE
@@ -35,14 +36,38 @@ namespace boost::math::tools {
 template<typename Real, typename Z = int64_t>
 class simple_continued_fraction {
 public:
-    simple_continued_fraction(Real x) : x_{x} {
+    typedef Z int_type;
+
+    simple_continued_fraction(std::vector<Z> 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;
         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.");
         }
+
+        if constexpr (std_precision == 2147483647) {
+           precision_ = x.backend().precision();
+        }
+
         b_.reserve(50);
         Real bj = floor(x);
         b_.push_back(static_cast<Z>(bj));
@@ -50,6 +75,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) {
@@ -57,12 +84,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<Real>::epsilon()*abs(x_))
-        {
+        const Real eps_abs_orig_x = std::numeric_limits<Real>::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<Z>(bj));
           x = 1/(x-bj);
@@ -77,16 +103,10 @@ 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_[b_.size() - 2] += 1;
-          b_.resize(b_.size() - 1);
-       }
-       b_.shrink_to_fit();
-       
-       for (size_t i = 1; i < b_.size(); ++i) {
+       canonicalize();
+
+       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 +121,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<Real>::quiet_NaN();
         }
          using std::log;
@@ -113,7 +134,7 @@ class simple_continued_fraction {
          // A random partial denominator has ~80% chance of being in this table:
          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) {
+         for (size_t i = 1; i < size_; ++i) {
             if (b_[i] < static_cast<Z>(logs.size())) {
                log_prod += logs[b_[i]];
             }
@@ -122,49 +143,57 @@ class simple_continued_fraction {
                log_prod += log(static_cast<Real>(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<Real>::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<Real>(b_[i]);
         }
         return n/denom;
     }
-    
+
+    // Note that this also includes the integer part (i.e. all the coefficients)
     const std::vector<Z>& partial_denominators() const {
       return b_;
     }
 
-    inline std::vector<Z>&& get_data() noexcept
-    {
+    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);
-
 private:
-    const Real x_;
+    static constexpr int std_precision = std::numeric_limits<Real>::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<Z> b_;
+
+    int precision_{std_precision};
 };
 
 
 template<typename Real, typename Z2>
 std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>& scf) {
-   constexpr const int p = std::numeric_limits<Real>::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)
    {
@@ -182,9 +211,38 @@ std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>&
 template<typename Real, typename Z = std::int64_t>
 inline auto simple_continued_fraction_coefficients(Real x)
 {
-    auto temp = simple_continued_fraction(x);
+    auto temp = simple_continued_fraction<Real, Z>(x);
     return temp.get_data();
 }
 
+// Can be used with `boost::rational` from <boost/rational.hpp>
+template <typename Rational, typename Real, typename Z = std::int64_t>
+inline Rational to_rational(const simple_continued_fraction<Real, Z>& 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<int_t>(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