From 4937338e159612c22b7c953d0c86eba4bef5b506 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:26:13 +0200 Subject: [PATCH] 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..c3eb52b891 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 (int 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