diff --git a/doc/equations/luroth_expansion.svg b/doc/equations/luroth_expansion.svg
new file mode 100644
index 000000000..1f7b4af5c
--- /dev/null
+++ b/doc/equations/luroth_expansion.svg
@@ -0,0 +1,2 @@
+
+
\ No newline at end of file
diff --git a/doc/internals/luroth_expansion.qbk b/doc/internals/luroth_expansion.qbk
new file mode 100644
index 000000000..70002dac4
--- /dev/null
+++ b/doc/internals/luroth_expansion.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:luroth_expansion Luroth Expansions]
+
+ #include
+ namespace boost::math::tools {
+
+ template
+ class luroth_expansion {
+ public:
+ luroth_expansion(Real x);
+
+ std::vector const & digits() const;
+
+ Real digit_geometric_mean() const;
+
+ template
+ friend std::ostream& operator<<(std::ostream& out, luroth_expansion& luroth);
+ };
+ }
+
+
+The `luroth_expansion` class provided by Boost expands a floating point number into a Lüroth representation, i.e.,
+
+[$../equations/luroth_expansion.svg]
+
+The numbers /d/[sub i] are called digits or denominators; we use the terminology digits, since technically in our notation /d/[sub 0] is not a denominator.
+
+Here's a minimal working example:
+
+ using boost::math::constants::pi;
+ using boost::math::tools::luroth_expansion;
+ auto luroth = luroth_expansion(pi());
+ std::cout << "π ≈ " << luroth << "\n";
+ // Prints:
+ // π ≈ ((3; 7, 1, 1, 1, 2, 1, 4, 23, 4, 1, 1, 1, 1, 80, 1, 1, 5))
+
+The class computes denominators while simultaneously computing convergents.
+Once a convergent is within a few ulps of the input value, the computation stops.
+
+/Nota bene:/ There is an alternative definition of the Lüroth representation where every digit is shifted by 1.
+We follow the definition given in Kalpazidou; with the modification that we do not constrain the input to be in the interval [0,1]
+and let the first digit be the floor of the input.
+
+For almost all real numbers, the geometric mean of the digits converges to a constant which is approximately 2.2001610580.
+This is "Khinchin's constant" for the Lüroth representation.
+
+
+
+[heading References]
+
+* Kalpazidou, Sofia. "Khintchine's constant for Lüroth representation." Journal of Number Theory 29.2 (1988): 196-205.
+
+* Finch, Steven R. Mathematical constants. Cambridge university press, 2003.
+
+
+[endsect]
diff --git a/doc/math.qbk b/doc/math.qbk
index df0770f6f..dd97e7228 100644
--- a/doc/math.qbk
+++ b/doc/math.qbk
@@ -753,6 +753,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22.
[include internals/fraction.qbk]
[include internals/simple_continued_fraction.qbk]
[include internals/centered_continued_fraction.qbk]
+[include internals/luroth_expansion.qbk]
[include internals/recurrence.qbk]
[/include internals/rational.qbk] [/moved to tools]
[include internals/tuple.qbk]
diff --git a/example/luroth.cpp b/example/luroth.cpp
new file mode 100644
index 000000000..cfaa55387
--- /dev/null
+++ b/example/luroth.cpp
@@ -0,0 +1,18 @@
+// (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
+#include
+#include
+
+using boost::math::constants::pi;
+using boost::math::tools::luroth_expansion;
+using boost::multiprecision::mpfr_float;
+
+int main() {
+ using Real = mpfr_float;
+ mpfr_float::default_precision(1024);
+ auto luroth = luroth_expansion(pi());
+ std::cout << luroth << "\n";
+}
diff --git a/include/boost/math/tools/luroth_expansion.hpp b/include/boost/math/tools/luroth_expansion.hpp
new file mode 100644
index 000000000..385a02e04
--- /dev/null
+++ b/include/boost/math/tools/luroth_expansion.hpp
@@ -0,0 +1,138 @@
+// (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_LUROTH_EXPANSION_HPP
+#define BOOST_MATH_TOOLS_LUROTH_EXPANSION_HPP
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace boost::math::tools {
+
+template
+class luroth_expansion {
+public:
+ luroth_expansion(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 a Lüroth representation.");
+ }
+ d_.reserve(50);
+ Real dn = floor(x);
+ d_.push_back(static_cast(dn));
+ if (dn == x)
+ {
+ d_.shrink_to_fit();
+ return;
+ }
+ // This attempts to follow the notation of:
+ // "Khinchine's constant for Luroth Representation", by Sophia Kalpazidou.
+ x = x - dn;
+ Real computed = dn;
+ Real prod = 1;
+ // Let the error bound grow by 1 ULP/iteration.
+ // I haven't done the error analysis to show that this is an expected rate of error growth,
+ // but if you don't do this, you can easily get into an infinite loop.
+ Real i = 1;
+ Real scale = std::numeric_limits::epsilon()*abs(x_)/2;
+ while (abs(x_ - computed) > (i++)*scale)
+ {
+ Real recip = 1/x;
+ Real dn = floor(recip);
+ // x = n + 1/k => lur(x) = ((n; k - 1))
+ // Note that this is a bit different than Kalpazidou (examine the half-open interval of definition carefully).
+ // One way to examine this definition is better for rationals (it never happens for irrationals)
+ // is to consider i + 1/3. If you follow Kalpazidou, then you get ((i, 3, 0)); a zero digit!
+ // That's bad since it destroys uniqueness and also breaks the computation of the geometric mean.
+ if (recip == dn) {
+ d_.push_back(static_cast(dn - 1));
+ break;
+ }
+ d_.push_back(static_cast(dn));
+ Real tmp = 1/(dn+1);
+ computed += prod*tmp;
+ prod *= tmp/dn;
+ x = dn*(dn+1)*(x - tmp);
+ }
+
+ for (size_t i = 1; i < d_.size(); ++i)
+ {
+ // Sanity check:
+ if (d_[i] <= 0)
+ {
+ throw std::domain_error("Found a digit <= 0; this is an error.");
+ }
+ }
+ d_.shrink_to_fit();
+ }
+
+
+ const std::vector& digits() const {
+ return d_;
+ }
+
+ // Under the assumption of 'randomness', this mean converges to 2.2001610580.
+ // See Finch, Mathematical Constants, section 1.8.1.
+ Real digit_geometric_mean() const {
+ if (d_.size() == 1) {
+ return std::numeric_limits::quiet_NaN();
+ }
+ using std::log;
+ using std::exp;
+ Real g = 0;
+ for (size_t i = 1; i < d_.size(); ++i) {
+ g += log(static_cast(d_[i]));
+ }
+ return exp(g/(d_.size() - 1));
+ }
+
+ template
+ friend std::ostream& operator<<(std::ostream& out, luroth_expansion& scf);
+
+private:
+ const Real x_;
+ std::vector d_;
+};
+
+
+template
+std::ostream& operator<<(std::ostream& out, luroth_expansion& luroth)
+{
+ constexpr const int p = std::numeric_limits::max_digits10;
+ if constexpr (p == 2147483647)
+ {
+ out << std::setprecision(luroth.x_.backend().precision());
+ }
+ else
+ {
+ out << std::setprecision(p);
+ }
+
+ out << "((" << luroth.d_.front();
+ if (luroth.d_.size() > 1)
+ {
+ out << "; ";
+ for (size_t i = 1; i < luroth.d_.size() -1; ++i)
+ {
+ out << luroth.d_[i] << ", ";
+ }
+ out << luroth.d_.back();
+ }
+ out << "))";
+ return out;
+}
+
+
+}
+#endif
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index c629dc290..a98f8dc31 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -896,6 +896,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" : -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" : -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" : -lquadmath ] ]
+ [ run luroth_expansion_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -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" : -lquadmath ] ]
diff --git a/test/luroth_expansion_test.cpp b/test/luroth_expansion_test.cpp
new file mode 100644
index 000000000..f9664f108
--- /dev/null
+++ b/test/luroth_expansion_test.cpp
@@ -0,0 +1,95 @@
+/*
+ * 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
+#include
+#ifdef BOOST_HAS_FLOAT128
+#include
+using boost::multiprecision::float128;
+#endif
+#include
+
+using boost::math::tools::luroth_expansion;
+using boost::multiprecision::cpp_bin_float_100;
+using boost::math::constants::pi;
+
+template
+void test_integral()
+{
+ for (int64_t i = -20; i < 20; ++i) {
+ Real ii = i;
+ auto luroth = luroth_expansion(ii);
+ auto const & a = luroth.digits();
+ CHECK_EQUAL(size_t(1), a.size());
+ CHECK_EQUAL(i, a.front());
+ }
+}
+
+
+template
+void test_halves()
+{
+ // x = n + 1/k => lur(x) = ((n; k - 1))
+ // Note that this is a bit different that Kalpazidou (examine the half-open interval of definition carefully).
+ // One way to examine this definition is correct for rationals (it never happens for irrationals)
+ // is to consider i + 1/3. If you follow Kalpazidou, then you get ((i, 3, 0)); a zero digit!
+ // That's bad since it destroys uniqueness and also breaks the computation of the geometric mean.
+ for (int64_t i = -20; i < 20; ++i) {
+ Real x = i + Real(1)/Real(2);
+ auto luroth = luroth_expansion(x);
+ auto const & a = luroth.digits();
+ CHECK_EQUAL(size_t(2), a.size());
+ CHECK_EQUAL(i, a.front());
+ CHECK_EQUAL(int64_t(1), a.back());
+ }
+
+ for (int64_t i = -20; i < 20; ++i) {
+ Real x = i + Real(1)/Real(4);
+ auto luroth = luroth_expansion(x);
+ auto const & a = luroth.digits();
+ CHECK_EQUAL(size_t(2), a.size());
+ CHECK_EQUAL(i, a.front());
+ CHECK_EQUAL(int64_t(3), a.back());
+ }
+
+ for (int64_t i = -20; i < 20; ++i) {
+ Real x = i + Real(1)/Real(8);
+ auto luroth = luroth_expansion(x);
+ auto const & a = luroth.digits();
+ CHECK_EQUAL(size_t(2), a.size());
+ CHECK_EQUAL(i, a.front());
+ CHECK_EQUAL(int64_t(7), a.back());
+ }
+ // 1/3 is a pain because it's not representable:
+ Real x = Real(1)/Real(3);
+ auto luroth = luroth_expansion(x);
+ auto const & a = luroth.digits();
+ CHECK_EQUAL(size_t(2), a.size());
+ CHECK_EQUAL(int64_t(0), a.front());
+ CHECK_EQUAL(int64_t(2), a.back());
+}
+
+
+int main()
+{
+ test_integral();
+ test_integral();
+ test_integral();
+ test_integral();
+
+ test_halves();
+ test_halves();
+ test_halves();
+ test_halves();
+
+ #ifdef BOOST_HAS_FLOAT128
+ test_integral();
+ test_halves();
+ #endif
+ return boost::math::test::report_errors();
+}