]>
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_CENTERED_CONTINUED_FRACTION_HPP | |
7 | #define BOOST_MATH_TOOLS_CENTERED_CONTINUED_FRACTION_HPP | |
8 | ||
1e59de90 TL |
9 | #include <cmath> |
10 | #include <cstdint> | |
20effc67 TL |
11 | #include <vector> |
12 | #include <ostream> | |
13 | #include <iomanip> | |
20effc67 TL |
14 | #include <limits> |
15 | #include <stdexcept> | |
1e59de90 TL |
16 | #include <sstream> |
17 | #include <array> | |
18 | #include <type_traits> | |
19 | #include <boost/math/tools/is_standalone.hpp> | |
20 | ||
21 | #ifndef BOOST_MATH_STANDALONE | |
20effc67 | 22 | #include <boost/core/demangle.hpp> |
1e59de90 | 23 | #endif |
20effc67 TL |
24 | |
25 | namespace boost::math::tools { | |
26 | ||
27 | template<typename Real, typename Z = int64_t> | |
28 | class centered_continued_fraction { | |
29 | public: | |
30 | centered_continued_fraction(Real x) : x_{x} { | |
31 | static_assert(std::is_integral_v<Z> && std::is_signed_v<Z>, | |
32 | "Centered continued fractions require signed integer types."); | |
33 | using std::round; | |
34 | using std::abs; | |
35 | using std::sqrt; | |
36 | using std::isfinite; | |
37 | if (!isfinite(x)) | |
38 | { | |
39 | throw std::domain_error("Cannot convert non-finites into continued fractions."); | |
40 | } | |
41 | b_.reserve(50); | |
42 | Real bj = round(x); | |
43 | b_.push_back(static_cast<Z>(bj)); | |
44 | if (bj == x) | |
45 | { | |
46 | b_.shrink_to_fit(); | |
47 | return; | |
48 | } | |
49 | x = 1/(x-bj); | |
50 | Real f = bj; | |
51 | if (bj == 0) | |
52 | { | |
1e59de90 | 53 | f = 16*(std::numeric_limits<Real>::min)(); |
20effc67 TL |
54 | } |
55 | Real C = f; | |
56 | Real D = 0; | |
57 | int i = 0; | |
58 | while (abs(f - x_) >= (1 + i++)*std::numeric_limits<Real>::epsilon()*abs(x_)) | |
59 | { | |
60 | bj = round(x); | |
61 | b_.push_back(static_cast<Z>(bj)); | |
62 | x = 1/(x-bj); | |
63 | D += bj; | |
64 | if (D == 0) { | |
1e59de90 | 65 | D = 16*(std::numeric_limits<Real>::min)(); |
20effc67 TL |
66 | } |
67 | C = bj + 1/C; | |
68 | if (C==0) | |
69 | { | |
1e59de90 | 70 | C = 16*(std::numeric_limits<Real>::min)(); |
20effc67 TL |
71 | } |
72 | D = 1/D; | |
73 | f *= (C*D); | |
74 | } | |
75 | // Deal with non-uniqueness of continued fractions: [a0; a1, ..., an, 1] = a0; a1, ..., an + 1]. | |
76 | if (b_.size() > 2 && b_.back() == 1) | |
77 | { | |
78 | b_[b_.size() - 2] += 1; | |
79 | b_.resize(b_.size() - 1); | |
80 | } | |
81 | b_.shrink_to_fit(); | |
82 | ||
83 | for (size_t i = 1; i < b_.size(); ++i) | |
84 | { | |
85 | if (b_[i] == 0) { | |
86 | std::ostringstream oss; | |
87 | oss << "Found a zero partial denominator: b[" << i << "] = " << b_[i] << "." | |
1e59de90 | 88 | #ifndef BOOST_MATH_STANDALONE |
20effc67 | 89 | << " This means the integer type '" << boost::core::demangle(typeid(Z).name()) |
1e59de90 TL |
90 | #else |
91 | << " This means the integer type '" << typeid(Z).name() | |
92 | #endif | |
20effc67 TL |
93 | << "' has overflowed and you need to use a wider type," |
94 | << " or there is a bug."; | |
95 | throw std::overflow_error(oss.str()); | |
96 | } | |
97 | } | |
98 | } | |
99 | ||
100 | Real khinchin_geometric_mean() const { | |
101 | if (b_.size() == 1) | |
102 | { | |
103 | return std::numeric_limits<Real>::quiet_NaN(); | |
104 | } | |
105 | using std::log; | |
106 | using std::exp; | |
107 | using std::abs; | |
108 | 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))}; | |
109 | Real log_prod = 0; | |
110 | for (size_t i = 1; i < b_.size(); ++i) | |
111 | { | |
112 | if (abs(b_[i]) < static_cast<Z>(logs.size())) | |
113 | { | |
114 | log_prod += logs[abs(b_[i])]; | |
115 | } | |
116 | else | |
117 | { | |
118 | log_prod += log(static_cast<Real>(abs(b_[i]))); | |
119 | } | |
120 | } | |
121 | log_prod /= (b_.size()-1); | |
122 | return exp(log_prod); | |
123 | } | |
124 | ||
125 | const std::vector<Z>& partial_denominators() const { | |
126 | return b_; | |
127 | } | |
128 | ||
129 | template<typename T, typename Z2> | |
130 | friend std::ostream& operator<<(std::ostream& out, centered_continued_fraction<T, Z2>& ccf); | |
131 | ||
132 | private: | |
133 | const Real x_; | |
134 | std::vector<Z> b_; | |
135 | }; | |
136 | ||
137 | ||
138 | template<typename Real, typename Z2> | |
139 | std::ostream& operator<<(std::ostream& out, centered_continued_fraction<Real, Z2>& scf) { | |
140 | constexpr const int p = std::numeric_limits<Real>::max_digits10; | |
141 | if constexpr (p == 2147483647) | |
142 | { | |
143 | out << std::setprecision(scf.x_.backend().precision()); | |
144 | } | |
145 | else | |
146 | { | |
147 | out << std::setprecision(p); | |
148 | } | |
149 | ||
150 | out << "[" << scf.b_.front(); | |
151 | if (scf.b_.size() > 1) | |
152 | { | |
153 | out << "; "; | |
154 | for (size_t i = 1; i < scf.b_.size() -1; ++i) | |
155 | { | |
156 | out << scf.b_[i] << ", "; | |
157 | } | |
158 | out << scf.b_.back(); | |
159 | } | |
160 | out << "]"; | |
161 | return out; | |
162 | } | |
163 | ||
164 | ||
165 | } | |
166 | #endif |