]>
Commit | Line | Data |
---|---|---|
20effc67 TL |
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_SIMPLE_CONTINUED_FRACTION_HPP | |
7 | #define BOOST_MATH_TOOLS_SIMPLE_CONTINUED_FRACTION_HPP | |
8 | ||
1e59de90 | 9 | #include <array> |
20effc67 TL |
10 | #include <vector> |
11 | #include <ostream> | |
12 | #include <iomanip> | |
13 | #include <cmath> | |
14 | #include <limits> | |
15 | #include <stdexcept> | |
1e59de90 TL |
16 | #include <sstream> |
17 | #include <boost/math/tools/is_standalone.hpp> | |
18 | ||
19 | #ifndef BOOST_MATH_STANDALONE | |
20effc67 | 20 | #include <boost/core/demangle.hpp> |
1e59de90 | 21 | #endif |
20effc67 TL |
22 | |
23 | namespace boost::math::tools { | |
24 | ||
25 | template<typename Real, typename Z = int64_t> | |
26 | class simple_continued_fraction { | |
27 | public: | |
28 | simple_continued_fraction(Real x) : x_{x} { | |
29 | using std::floor; | |
30 | using std::abs; | |
31 | using std::sqrt; | |
32 | using std::isfinite; | |
33 | if (!isfinite(x)) { | |
34 | throw std::domain_error("Cannot convert non-finites into continued fractions."); | |
35 | } | |
36 | b_.reserve(50); | |
37 | Real bj = floor(x); | |
38 | b_.push_back(static_cast<Z>(bj)); | |
39 | if (bj == x) { | |
40 | b_.shrink_to_fit(); | |
41 | return; | |
42 | } | |
43 | x = 1/(x-bj); | |
44 | Real f = bj; | |
45 | if (bj == 0) { | |
1e59de90 | 46 | f = 16*(std::numeric_limits<Real>::min)(); |
20effc67 TL |
47 | } |
48 | Real C = f; | |
49 | Real D = 0; | |
50 | int i = 0; | |
51 | // the "1 + i++" lets the error bound grow slowly with the number of convergents. | |
52 | // I have not worked out the error propagation of the Modified Lentz's method to see if it does indeed grow at this rate. | |
53 | // Numerical Recipes claims that no one has worked out the error analysis of the modified Lentz's method. | |
54 | while (abs(f - x_) >= (1 + i++)*std::numeric_limits<Real>::epsilon()*abs(x_)) | |
55 | { | |
56 | bj = floor(x); | |
57 | b_.push_back(static_cast<Z>(bj)); | |
58 | x = 1/(x-bj); | |
59 | D += bj; | |
60 | if (D == 0) { | |
1e59de90 | 61 | D = 16*(std::numeric_limits<Real>::min)(); |
20effc67 TL |
62 | } |
63 | C = bj + 1/C; | |
64 | if (C==0) { | |
1e59de90 | 65 | C = 16*(std::numeric_limits<Real>::min)(); |
20effc67 TL |
66 | } |
67 | D = 1/D; | |
68 | f *= (C*D); | |
69 | } | |
70 | // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1]. | |
71 | // The shorter representation is considered the canonical representation, | |
72 | // so if we compute a non-canonical representation, change it to canonical: | |
73 | if (b_.size() > 2 && b_.back() == 1) { | |
74 | b_[b_.size() - 2] += 1; | |
75 | b_.resize(b_.size() - 1); | |
76 | } | |
77 | b_.shrink_to_fit(); | |
78 | ||
79 | for (size_t i = 1; i < b_.size(); ++i) { | |
80 | if (b_[i] <= 0) { | |
81 | std::ostringstream oss; | |
82 | oss << "Found a negative partial denominator: b[" << i << "] = " << b_[i] << "." | |
1e59de90 | 83 | #ifndef BOOST_MATH_STANDALONE |
20effc67 | 84 | << " This means the integer type '" << boost::core::demangle(typeid(Z).name()) |
1e59de90 TL |
85 | #else |
86 | << " This means the integer type '" << typeid(Z).name() | |
87 | #endif | |
20effc67 TL |
88 | << "' has overflowed and you need to use a wider type," |
89 | << " or there is a bug."; | |
90 | throw std::overflow_error(oss.str()); | |
91 | } | |
92 | } | |
93 | } | |
94 | ||
95 | Real khinchin_geometric_mean() const { | |
96 | if (b_.size() == 1) { | |
97 | return std::numeric_limits<Real>::quiet_NaN(); | |
98 | } | |
99 | using std::log; | |
100 | using std::exp; | |
101 | // Precompute the most probable logarithms. See the Gauss-Kuzmin distribution for details. | |
1e59de90 | 102 | // Example: b_i = 1 has probability -log_2(3/4) ~ .415: |
20effc67 TL |
103 | // A random partial denominator has ~80% chance of being in this table: |
104 | 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))}; | |
105 | Real log_prod = 0; | |
106 | for (size_t i = 1; i < b_.size(); ++i) { | |
107 | if (b_[i] < static_cast<Z>(logs.size())) { | |
108 | log_prod += logs[b_[i]]; | |
109 | } | |
110 | else | |
111 | { | |
112 | log_prod += log(static_cast<Real>(b_[i])); | |
113 | } | |
114 | } | |
115 | log_prod /= (b_.size()-1); | |
116 | return exp(log_prod); | |
117 | } | |
118 | ||
119 | Real khinchin_harmonic_mean() const { | |
120 | if (b_.size() == 1) { | |
121 | return std::numeric_limits<Real>::quiet_NaN(); | |
122 | } | |
123 | Real n = b_.size() - 1; | |
124 | Real denom = 0; | |
125 | for (size_t i = 1; i < b_.size(); ++i) { | |
126 | denom += 1/static_cast<Real>(b_[i]); | |
127 | } | |
128 | return n/denom; | |
129 | } | |
130 | ||
131 | const std::vector<Z>& partial_denominators() const { | |
132 | return b_; | |
133 | } | |
134 | ||
135 | template<typename T, typename Z2> | |
136 | friend std::ostream& operator<<(std::ostream& out, simple_continued_fraction<T, Z2>& scf); | |
137 | ||
138 | private: | |
139 | const Real x_; | |
140 | std::vector<Z> b_; | |
141 | }; | |
142 | ||
143 | ||
144 | template<typename Real, typename Z2> | |
145 | std::ostream& operator<<(std::ostream& out, simple_continued_fraction<Real, Z2>& scf) { | |
146 | constexpr const int p = std::numeric_limits<Real>::max_digits10; | |
147 | if constexpr (p == 2147483647) { | |
148 | out << std::setprecision(scf.x_.backend().precision()); | |
149 | } else { | |
150 | out << std::setprecision(p); | |
151 | } | |
152 | ||
153 | out << "[" << scf.b_.front(); | |
154 | if (scf.b_.size() > 1) | |
155 | { | |
156 | out << "; "; | |
157 | for (size_t i = 1; i < scf.b_.size() -1; ++i) | |
158 | { | |
159 | out << scf.b_[i] << ", "; | |
160 | } | |
161 | out << scf.b_.back(); | |
162 | } | |
163 | out << "]"; | |
164 | return out; | |
165 | } | |
166 | ||
167 | ||
168 | } | |
169 | #endif |