Centered continued fractions (#379)

* Centered continued fraction [CI SKIP]

* Document centered cfrac. [CI SKIP]

* Unit tests for centered continued fraction [CI SKIP]

* Kick off build.

* Fix syntax error in docs [CI SKIP]

* Fix ADL.

doc/internals/centered_continued_fraction.qbk

@@ -0,0 +1,55 @@
+[/
+  Copyright Nick Thompson, 2020
+  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).
+]
+
+[section:centered_continued_fraction Centered Continued Fractions]
+
+    #include <boost/math/tools/centered_continued_fraction.hpp>
+    namespace boost::math::tools {
+
+    template<typename Real>
+    class centered_continued_fraction {
+    public:
+        centered_continued_fraction(Real x);
+
+        Real khinchin_geometric_mean() const;
+
+        template<typename T, typename Z_>
+        friend std::ostream& operator<<(std::ostream& out, centered_continued_fraction<T, Z>& scf);
+    };
+    }
+
+
+The `centered_continued_fraction` class provided by Boost affords the ability to convert a floating point number into a centered continued fraction.
+Unlike the simple continued fraction, a centered continued fraction allows partial denominators to be negative.
+
+Here's a minimal working example:
+
+    using boost::math::constants::pi;
+    using boost::math::tools::centered_continued_fraction;
+    auto cfrac = centered_continued_fraction(pi<long double>());
+    std::cout << "π ≈ " << cfrac << "\n";
+    // Prints:
+    // π ≈ [3; 7, 16, -294, 3, -4, 5, -15, -3, 2, 2]
+
+This function achieves the same end as the Maple syntax
+
+[$../images/maple_cfrac.png]
+
+You should convince yourself that the Maple output and the notation we use are in fact the same thing.
+
+The class computes partial denominators which simultaneously computing convergents with the modified Lentz's algorithm.
+Once a convergent is within a few ulps of the input value, the computation stops.
+
+Just like simple continued fractions, centered continued fractions admit a "Khinchin constant".
+The value of this constant is ~5.454517 (see [@http://jeremiebourdon.free.fr/data/Khintchine.pdf here]).
+We can evaluate the "empirical Khinchin constant" for a particular number via
+
+    cfrac.khinchin_geometric_mean();
+
+If this is close to 5.454517, then the number is probably uninteresting with respect to its characterization as a centered continued fraction.
+
+[endsect]

doc/math.qbk

@@ -752,6 +752,7 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include internals/agm.qbk]
 [include internals/fraction.qbk]
 [include internals/simple_continued_fraction.qbk]
+[include internals/centered_continued_fraction.qbk]
 [include internals/recurrence.qbk]
 [/include internals/rational.qbk] [/moved to tools]
 [include internals/tuple.qbk]

example/centered_continued_fraction.cpp

@@ -0,0 +1,46 @@
+//  (C) Copyright Nick Thompson 2020.
+//  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)
+#include <iostream>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/tools/centered_continued_fraction.hpp>
+#include <boost/multiprecision/mpfr.hpp>
+
+using boost::math::constants::root_two;
+using boost::math::constants::phi;
+using boost::math::constants::pi;
+using boost::math::constants::e;
+using boost::math::constants::zeta_three;
+using boost::math::tools::centered_continued_fraction;
+using boost::multiprecision::mpfr_float;
+
+int main()
+{
+    using Real = mpfr_float;
+    int p = 10000;
+    mpfr_float::default_precision(p);
+    auto phi_cfrac = centered_continued_fraction(phi<Real>());
+    std::cout << "φ ≈ " << phi_cfrac << "\n";
+    std::cout << "Khinchin mean: " << std::setprecision(10) << phi_cfrac.khinchin_geometric_mean() << "\n\n\n";
+
+    auto pi_cfrac = centered_continued_fraction(pi<Real>());
+    std::cout << "π ≈ " << pi_cfrac << "\n";
+    std::cout << "Khinchin mean: " << std::setprecision(10) << pi_cfrac.khinchin_geometric_mean() << "\n\n\n";
+
+    auto rt_cfrac = centered_continued_fraction(root_two<Real>());
+    std::cout << "√2 ≈ " << rt_cfrac << "\n";
+    std::cout << "Khinchin mean: " << std::setprecision(10) <<  rt_cfrac.khinchin_geometric_mean() << "\n\n\n";
+
+    auto e_cfrac = centered_continued_fraction(e<Real>());
+    std::cout << "e ≈ " << e_cfrac << "\n";
+    std::cout << "Khinchin mean: " << std::setprecision(10) <<  e_cfrac.khinchin_geometric_mean() << "\n\n\n";
+
+    auto z_cfrac = centered_continued_fraction(zeta_three<Real>());
+    std::cout << "ζ(3) ≈ " << z_cfrac << "\n";
+    std::cout << "Khinchin mean: " << std::setprecision(10) <<  z_cfrac.khinchin_geometric_mean() << "\n\n\n";
+
+
+    // http://jeremiebourdon.free.fr/data/Khintchine.pdf
+    std::cout << "The expected Khinchin mean for a random centered continued fraction is 5.45451724454\n";
+}

include/boost/math/tools/centered_continued_fraction.hpp

@@ -0,0 +1,154 @@
+//  (C) Copyright Nick Thompson 2020.
+//  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)
+
+#ifndef BOOST_MATH_TOOLS_CENTERED_CONTINUED_FRACTION_HPP
+#define BOOST_MATH_TOOLS_CENTERED_CONTINUED_FRACTION_HPP
+
+#include <vector>
+#include <ostream>
+#include <iomanip>
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+#include <boost/core/demangle.hpp>
+
+namespace boost::math::tools {
+
+template<typename Real, typename Z = int64_t>
+class centered_continued_fraction {
+public:
+    centered_continued_fraction(Real x) : x_{x} {
+        static_assert(std::is_integral_v<Z> && std::is_signed_v<Z>,
+                      "Centered continued fractions require signed integer types.");
+        using std::round;
+        using std::abs;
+        using std::sqrt;
+        using std::isfinite;
+        if (!isfinite(x))
+        {
+            throw std::domain_error("Cannot convert non-finites into continued fractions.");  
+        }
+        b_.reserve(50);
+        Real bj = round(x);
+        b_.push_back(static_cast<Z>(bj));
+        if (bj == x)
+        {
+            b_.shrink_to_fit();
+            return;
+        }
+        x = 1/(x-bj);
+        Real f = bj;
+        if (bj == 0)
+        {
+            f = 16*std::numeric_limits<Real>::min();
+        }
+        Real C = f;
+        Real D = 0;
+        int i = 0;
+        while (abs(f - x_) >= (1 + i++)*std::numeric_limits<Real>::epsilon()*abs(x_))
+        {
+            bj = round(x);
+            b_.push_back(static_cast<Z>(bj));
+            x = 1/(x-bj);
+            D += bj;
+            if (D == 0) {
+                D = 16*std::numeric_limits<Real>::min();
+            }
+            C = bj + 1/C;
+            if (C==0)
+            {
+                C = 16*std::numeric_limits<Real>::min();
+            }
+            D = 1/D;
+            f *= (C*D);
+        }
+        // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1].
+        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)
+        {
+            if (b_[i] == 0) {
+                std::ostringstream oss;
+                oss << "Found a zero partial denominator: b[" << i << "] = " << b_[i] << "."
+                    << " This means the integer type '" << boost::core::demangle(typeid(Z).name())
+                    << "' has overflowed and you need to use a wider type,"
+                    << " or there is a bug.";
+                throw std::overflow_error(oss.str());
+            }
+        }
+    }
+
+    Real khinchin_geometric_mean() const {
+        if (b_.size() == 1)
+        { 
+            return std::numeric_limits<Real>::quiet_NaN();
+        }
+        using std::log;
+        using std::exp;
+        using std::abs;
+        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))};
+        Real log_prod = 0;
+        for (size_t i = 1; i < b_.size(); ++i)
+        {
+            if (abs(b_[i]) < static_cast<Z>(logs.size()))
+            {
+                log_prod += logs[abs(b_[i])];
+            }
+            else
+            {
+                log_prod += log(static_cast<Real>(abs(b_[i])));
+            }
+        }
+        log_prod /= (b_.size()-1);
+        return exp(log_prod);
+    }
+
+    const std::vector<Z>& partial_denominators() const {
+        return b_;
+    }
+    
+    template<typename T, typename Z2>
+    friend std::ostream& operator<<(std::ostream& out, centered_continued_fraction<T, Z2>& ccf);
+
+private:
+    const Real x_;
+    std::vector<Z> b_;
+};
+
+
+template<typename Real, typename Z2>
+std::ostream& operator<<(std::ostream& out, centered_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 << "[" << scf.b_.front();
+    if (scf.b_.size() > 1)
+    {
+        out << "; ";
+        for (size_t i = 1; i < scf.b_.size() -1; ++i)
+        {
+            out << scf.b_[i] << ", ";
+        }
+        out << scf.b_.back();
+    }
+    out << "]";
+    return out;
+}
+
+
+}
+#endif

include/boost/math/tools/simple_continued_fraction.hpp

@@ -140,7 +140,7 @@ std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>&
       out << std::setprecision(p);
    }
    
-   out << scf.x_ << " ≈ [" << scf.b_.front();
+   out << "[" << scf.b_.front();
    if (scf.b_.size() > 1)
    {
       out << "; ";

test/Jamfile.v2

@@ -895,6 +895,7 @@ test-suite misc :
     ]
    [ run test_constants.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
    [ run simple_continued_fraction_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
+   [ run centered_continued_fraction_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
    [ run test_classify.cpp pch ../../test/build//boost_unit_test_framework  ]
    [ run test_error_handling.cpp ../../test/build//boost_unit_test_framework  ]
    [ run legendre_stieltjes_test.cpp ../../test/build//boost_unit_test_framework : : : [ requires cxx11_auto_declarations cxx11_range_based_for ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ]  ]

test/centered_continued_fraction_test.cpp

@@ -0,0 +1,148 @@
+/*
+ * Copyright Nick Thompson, 2020
+ * 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)
+ */
+
+#include "math_unit_test.hpp"
+#include <boost/math/tools/centered_continued_fraction.hpp>
+#include <boost/math/constants/constants.hpp>
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+using boost::multiprecision::float128;
+#endif
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+using boost::math::tools::centered_continued_fraction;
+using boost::multiprecision::cpp_bin_float_100;
+using boost::math::constants::pi;
+
+template<class Real>
+void test_integral()
+{
+    for (int64_t i = -20; i < 20; ++i) {
+        Real ii = i;
+        auto cfrac = centered_continued_fraction<Real>(ii);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(1), a.size());
+        CHECK_EQUAL(i, a.front());
+    }
+}
+
+template<class Real>
+void test_halves()
+{
+    for (int64_t i = 0; i < 20; ++i) {
+        // std::round rounds 0.5 up if i > 0, and rounds down if i < 0.
+        // Should we care? These representations are not unique anyway;
+        // In any case, this behavior agrees with Maple.
+        Real x = i + Real(1)/Real(2);
+        auto cfrac = centered_continued_fraction<Real>(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(2), a.size());
+        CHECK_EQUAL(i + 1, a.front());
+        CHECK_EQUAL(int64_t(-2), a.back());
+    }
+
+    // We'll also test quarters; why not?
+    for (int64_t i = -20; i < 20; ++i) {
+        Real x = i + Real(1)/Real(4);
+        auto cfrac = centered_continued_fraction<Real>(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(2), a.size());
+        CHECK_EQUAL(i, a.front());
+        CHECK_EQUAL(int64_t(4), a.back());
+    }
+
+    for (int64_t i = -20; i < 20; ++i) {
+        Real x = i + Real(1)/Real(8);
+        auto cfrac = centered_continued_fraction<Real>(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(2), a.size());
+        CHECK_EQUAL(i, a.front());
+        CHECK_EQUAL(int64_t(8), a.back());
+    }
+
+    for (int64_t i = -20; i < 20; ++i) {
+        Real x = i + Real(3)/Real(4);
+        auto cfrac = centered_continued_fraction<Real>(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(2), a.size());
+        CHECK_EQUAL(i+1, a.front());
+        CHECK_EQUAL(int64_t(-4), a[1]);
+    }
+
+    for (int64_t i = -20; i < 20; ++i) {
+        Real x = i + Real(7)/Real(8);
+        auto cfrac = centered_continued_fraction<Real>(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(2), a.size());
+        CHECK_EQUAL(i + 1, a.front());
+        CHECK_EQUAL(int64_t(-8), a[1]);
+    }
+}
+
+template<typename Real>
+void test_simple()
+{
+    std::cout << "Testing rational numbers on type " << boost::core::demangle(typeid(Real).name()) << "\n";
+    {
+        Real x = Real(649)/200;
+        // ContinuedFraction[649/200] = [3; 4, 12, 4]
+        auto cfrac = centered_continued_fraction(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(4), a.size());
+        CHECK_EQUAL(int64_t(3), a[0]);
+        CHECK_EQUAL(int64_t(4), a[1]);
+        CHECK_EQUAL(int64_t(12), a[2]);
+        CHECK_EQUAL(int64_t(4), a[3]);
+    }
+
+    {
+        Real x = Real(415)/Real(93);
+        // [4; 2, 6, 7]:
+        auto cfrac = centered_continued_fraction(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(4), a.size());
+        CHECK_EQUAL(int64_t(4), a[0]);
+        CHECK_EQUAL(int64_t(2), a[1]);
+        CHECK_EQUAL(int64_t(6), a[2]);
+        CHECK_EQUAL(int64_t(7), a[3]);
+    }
+
+}
+
+template<typename Real>
+void test_khinchin()
+{
+    using std::sqrt;
+    auto rt_cfrac = centered_continued_fraction(sqrt(static_cast<Real>(2)));
+    Real K0 = rt_cfrac.khinchin_geometric_mean();
+    CHECK_ULP_CLOSE(Real(2), K0, 10);
+}
+
+
+int main()
+{
+    test_integral<float>();
+    test_integral<double>();
+    test_integral<long double>();
+    test_integral<cpp_bin_float_100>();
+    test_halves<float>();
+    test_halves<double>();
+    test_halves<long double>();
+    test_halves<cpp_bin_float_100>();
+    test_simple<float>();
+    test_simple<double>();
+    test_simple<long double>();
+    test_simple<cpp_bin_float_100>();
+    test_khinchin<cpp_bin_float_100>();
+    #ifdef BOOST_HAS_FLOAT128
+    test_integral<float128>();
+    test_halves<float128>();
+    test_simple<float128>();
+    test_khinchin<float128>();
+    #endif
+    return boost::math::test::report_errors();
+}