Simple continued fraction (#377)

* Simple continued fraction [CI SKIP]

* Comments on error analysis [CI SKIP]

* Simple continued fraction [CI SKIP]

* Clarify comment and kick off build.

doc/internals/simple_continued_fraction.qbk

@@ -0,0 +1,62 @@
+[/
+  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:simple_continued_fraction Simple Continued Fractions]
+
+    #include <boost/math/tools/simple_continued_fraction.hpp>
+    namespace boost::math::tools {
+
+    template<typename Real>
+    class simple_continued_fraction {
+    public:
+        simple_continued_fraction(Real x);
+        
+        Real khinchin_geometric_mean() const;
+        
+        Real khinchin_harmonic_mean() const;
+
+        template<typename T, typename Z_>
+        friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction<T, Z>& scf);
+    };
+    }
+
+
+The `simple_continued_fraction` class provided by Boost affords the ability to convert a floating point number into a simple continued fraction.
+In addition, we can answer a few questions about the number in question using this representation.
+
+Here's a minimal working example:
+
+    using boost::math::constants::pi;
+    using boost::math::tools::simple_continued_fraction;
+    auto cfrac = simple_continued_fraction(pi<long double>());
+    std::cout << "π ≈ " << cfrac << "\n";
+    // Prints:
+    // π ≈ 3.14159265358979323851 ≈ [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2]
+
+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.
+
+Note that every floating point number is a rational number, and this exact rational can be exactly converted to a finite continued fraction.
+This is perfectly sensible behavior, but we do not do it here.
+This is because when examining known values like π, it creates a large number of incorrect partial denominators, even if every bit of the binary representation is correct.
+
+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
+
+[$../equations/khinchin_geometric.svg]
+
+and Khinchin harmonic mean
+
+[$../equations/khinchin_harmonic.svg]
+
+If these approach Khinchin's constant /K/[sub 0] and /K/[sub -1] as the number of partial denominators goes to infinity, then our number is "uninteresting" with respect to the characterization.
+These violations are washed out if too many incorrect partial denominators are included in the expansion.
+
+Note: The convergence of these means to the Khinchin limit is exceedingly slow; we've used 30,000 decimal digits of π and only found two digits of agreement with /K/[sub 0].
+However, clear violations of are obvious, such as the continued fraction expansion of √2, whose Khinchin geometric mean is precisely 2.
+
+[endsect]

doc/math.qbk

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

example/to_continued_fraction.cpp

@@ -0,0 +1,36 @@
+//  (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/simple_continued_fraction.hpp>
+#include <boost/multiprecision/cpp_bin_float.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::simple_continued_fraction;
+
+int main()
+{
+    using Real = boost::multiprecision::cpp_bin_float_100;
+    auto phi_cfrac = simple_continued_fraction(phi<Real>());
+    std::cout << "φ ≈ " << phi_cfrac << "\n\n";
+
+    auto pi_cfrac = simple_continued_fraction(pi<Real>());
+    std::cout << "π ≈ " << pi_cfrac << "\n";
+    std::cout << "Known: [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 1, 1, 15, 3, 13, 1, 4, 2, 6, 6, 99, 1, 2, 2, 6, 3, 5, 1, 1, 6, 8, 1, 7, 1, 2, 3, 7, 1, 2, 1, 1, 12, 1, 1, 1, 3, 1, 1, 8, 1, 1, 2, 1, 6, 1, 1, 5, 2, 2, 3, 1, 2, 4, 4, 16, 1, 161, 45, 1, 22, 1, 2, 2, 1, 4, 1, 2, 24, 1, 2, 1, 3, 1, 2, 1, ...]\n\n";
+
+    auto rt_cfrac = simple_continued_fraction(root_two<Real>());
+    std::cout << "√2 ≈ " << rt_cfrac << "\n\n";
+
+    auto e_cfrac = simple_continued_fraction(e<Real>());
+    std::cout << "e ≈ " << e_cfrac << "\n";
+
+    // Correctness can be checked in Mathematica via: ContinuedFraction[Zeta[3], 500]
+    auto z_cfrac = simple_continued_fraction(zeta_three<Real>());
+    std::cout << "ζ(3) ≈ " << z_cfrac << "\n";
+}

include/boost/math/tools/simple_continued_fraction.hpp

@@ -0,0 +1,159 @@
+//  (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_SIMPLE_CONTINUED_FRACTION_HPP
+#define BOOST_MATH_TOOLS_SIMPLE_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 simple_continued_fraction {
+public:
+    simple_continued_fraction(Real x) : x_{x} {
+        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.");  
+        }
+        b_.reserve(50);
+        Real bj = floor(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;
+        // 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_))
+        {
+          bj = floor(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].
+       // 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) {
+         if (b_[i] <= 0) {
+            std::ostringstream oss;
+            oss << "Found a negative 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;
+         // 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))};
+         Real log_prod = 0;
+         for (size_t i = 1; i < b_.size(); ++i) {
+            if (b_[i] < static_cast<Z>(logs.size())) {
+               log_prod += logs[b_[i]];
+            }
+            else
+            {
+               log_prod += log(static_cast<Real>(b_[i]));
+            }
+         }
+         log_prod /= (b_.size()-1);
+         return exp(log_prod);
+    }
+    
+    Real khinchin_harmonic_mean() const {
+        if (b_.size() == 1) {
+          return std::numeric_limits<Real>::quiet_NaN();
+        }
+        Real n = b_.size() - 1;
+        Real denom = 0;
+        for (size_t i = 1; i < b_.size(); ++i) {
+            denom += 1/static_cast<Real>(b_[i]);
+        }
+        return n/denom;
+    }
+    
+    const std::vector<Z>& partial_denominators() const {
+      return b_;
+    }
+    
+    template<typename T, typename Z2>
+    friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction<T, Z2>& scf);
+
+private:
+    const Real x_;
+    std::vector<Z> b_;
+};
+
+
+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 << scf.x_ << " ≈ [" << 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 << "]\n";
+   return out;
+}
+
+
+}
+#endif

test/Jamfile.v2

@@ -894,6 +894,7 @@ test-suite misc :
       test_tr1_c_long_double
     ]
    [ 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 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/math_unit_test.hpp

@@ -323,9 +323,15 @@ bool check_equal(Real x, Real y, std::string const & filename, std::string const
   if (x != y) {
     std::ios_base::fmtflags f( std::cerr.flags() );
     std::cerr << "\033[0;31mError at " << filename << ":" << function << ":" << line << "\n";
-    std::cerr << "\033[0m  Expected " << x << " == " << y << " but:\n";
-    std::cerr << "  Expected =  " << std::defaultfloat << std::fixed << x << " = " << std::scientific << x << std::hexfloat << " = " << x << "\n";
-    std::cerr << "  Computed =  " << std::defaultfloat << std::fixed << y << " = " << std::scientific << y << std::hexfloat << " = " << y << "\n";
+    std::cerr << "\033[0m  Condition '" << x << " == " << y << "' is not satisfied:\n";
+    if (std::is_floating_point<Real>::value) {
+      std::cerr << "  Expected =  " << std::defaultfloat << std::fixed << x << " = " << std::scientific << x << std::hexfloat << " = " << x << "\n";
+      std::cerr << "  Computed =  " << std::defaultfloat << std::fixed << y << " = " << std::scientific << y << std::hexfloat << " = " << y << "\n";
+    } else {
+      std::cerr << "  Expected: " << x << " = " << "0x" << std::hex << x << "\n";
+      std::cerr << std::dec;
+      std::cerr << "  Computed: " << y << " = " << "0x" << std::hex << y << "\n";
+    }
     std::cerr.flags(f);
     ++detail::global_error_count;
     return false;

test/simple_continued_fraction_test.cpp

@@ -0,0 +1,160 @@
+/*
+ * 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/simple_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::simple_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 = simple_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 = -20; i < 20; ++i) {
+        Real x = i + Real(1)/Real(2);
+        auto cfrac = simple_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(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 = simple_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 = simple_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 = simple_continued_fraction<Real>(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(3), a.size());
+        CHECK_EQUAL(i, a.front());
+        CHECK_EQUAL(int64_t(1), a[1]);
+        CHECK_EQUAL(int64_t(3), a.back());
+    }
+
+    for (int64_t i = -20; i < 20; ++i) {
+        Real x = i + Real(7)/Real(8);
+        auto cfrac = simple_continued_fraction<Real>(x);
+        auto const & a = cfrac.partial_denominators();
+        CHECK_EQUAL(size_t(3), a.size());
+        CHECK_EQUAL(i, a.front());
+        CHECK_EQUAL(int64_t(1), a[1]);
+        CHECK_EQUAL(int64_t(7), a.back());
+    }
+}
+
+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 = simple_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 = simple_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()
+{
+    // These are simply sanity checks; the convergence is too slow otherwise:
+    auto cfrac = simple_continued_fraction(pi<Real>());
+    auto K0 = cfrac.khinchin_geometric_mean();
+    CHECK_MOLLIFIED_CLOSE(Real(2.6854520010), K0, 0.1);
+    auto Km1 = cfrac.khinchin_harmonic_mean();
+    CHECK_MOLLIFIED_CLOSE(Real(1.74540566240), Km1, 0.1);
+    
+    using std::sqrt;
+    auto rt_cfrac = simple_continued_fraction(sqrt(static_cast<Real>(2)));
+    K0 = rt_cfrac.khinchin_geometric_mean();
+    CHECK_ULP_CLOSE(Real(2), K0, 10);
+    Km1 = rt_cfrac.khinchin_harmonic_mean();
+    CHECK_ULP_CLOSE(Real(2), Km1, 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();
+}