1 // (C) Copyright Nick Thompson 2020.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #ifndef BOOST_MATH_TOOLS_LUROTH_EXPANSION_HPP
7 #define BOOST_MATH_TOOLS_LUROTH_EXPANSION_HPP
8
9 #include <vector>
10 #include <ostream>
11 #include <iomanip>
12 #include <cmath>
13 #include <limits>
14 #include <stdexcept>
15
16 namespace boost::math::tools {
17
18 template<typename Real, typename Z = int64_t>
19 class luroth_expansion {
20 public:
luroth_expansion(Real x)21 luroth_expansion(Real x) : x_{x}
22 {
23 using std::floor;
24 using std::abs;
25 using std::sqrt;
26 using std::isfinite;
27 if (!isfinite(x))
28 {
29 throw std::domain_error("Cannot convert non-finites into a Lüroth representation.");
30 }
31 d_.reserve(50);
32 Real dn = floor(x);
33 d_.push_back(static_cast<Z>(dn));
34 if (dn == x)
35 {
36 d_.shrink_to_fit();
37 return;
38 }
39 // This attempts to follow the notation of:
40 // "Khinchine's constant for Luroth Representation", by Sophia Kalpazidou.
41 x = x - dn;
42 Real computed = dn;
43 Real prod = 1;
44 // Let the error bound grow by 1 ULP/iteration.
45 // I haven't done the error analysis to show that this is an expected rate of error growth,
46 // but if you don't do this, you can easily get into an infinite loop.
47 Real i = 1;
48 Real scale = std::numeric_limits<Real>::epsilon()*abs(x_)/2;
49 while (abs(x_ - computed) > (i++)*scale)
50 {
51 Real recip = 1/x;
52 Real dn = floor(recip);
53 // x = n + 1/k => lur(x) = ((n; k - 1))
54 // Note that this is a bit different than Kalpazidou (examine the half-open interval of definition carefully).
55 // One way to examine this definition is better for rationals (it never happens for irrationals)
56 // is to consider i + 1/3. If you follow Kalpazidou, then you get ((i, 3, 0)); a zero digit!
57 // That's bad since it destroys uniqueness and also breaks the computation of the geometric mean.
58 if (recip == dn) {
59 d_.push_back(static_cast<Z>(dn - 1));
60 break;
61 }
62 d_.push_back(static_cast<Z>(dn));
63 Real tmp = 1/(dn+1);
64 computed += prod*tmp;
65 prod *= tmp/dn;
66 x = dn*(dn+1)*(x - tmp);
67 }
68
69 for (size_t i = 1; i < d_.size(); ++i)
70 {
71 // Sanity check:
72 if (d_[i] <= 0)
73 {
74 throw std::domain_error("Found a digit <= 0; this is an error.");
75 }
76 }
77 d_.shrink_to_fit();
78 }
79
80
digits() const81 const std::vector<Z>& digits() const {
82 return d_;
83 }
84
85 // Under the assumption of 'randomness', this mean converges to 2.2001610580.
86 // See Finch, Mathematical Constants, section 1.8.1.
digit_geometric_mean() const87 Real digit_geometric_mean() const {
88 if (d_.size() == 1) {
89 return std::numeric_limits<Real>::quiet_NaN();
90 }
91 using std::log;
92 using std::exp;
93 Real g = 0;
94 for (size_t i = 1; i < d_.size(); ++i) {
95 g += log(static_cast<Real>(d_[i]));
96 }
97 return exp(g/(d_.size() - 1));
98 }
99
100 template<typename T, typename Z2>
101 friend std::ostream& operator<<(std::ostream& out, luroth_expansion<T, Z2>& scf);
102
103 private:
104 const Real x_;
105 std::vector<Z> d_;
106 };
107
108
109 template<typename Real, typename Z2>
operator <<(std::ostream & out,luroth_expansion<Real,Z2> & luroth)110 std::ostream& operator<<(std::ostream& out, luroth_expansion<Real, Z2>& luroth)
111 {
112 constexpr const int p = std::numeric_limits<Real>::max_digits10;
113 if constexpr (p == 2147483647)
114 {
115 out << std::setprecision(luroth.x_.backend().precision());
116 }
117 else
118 {
119 out << std::setprecision(p);
120 }
121
122 out << "((" << luroth.d_.front();
123 if (luroth.d_.size() > 1)
124 {
125 out << "; ";
126 for (size_t i = 1; i < luroth.d_.size() -1; ++i)
127 {
128 out << luroth.d_[i] << ", ";
129 }
130 out << luroth.d_.back();
131 }
132 out << "))";
133 return out;
134 }
135
136
137 }
138 #endif
139