]>
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_ENGEL_EXPANSION_HPP | |
7 | #define BOOST_MATH_TOOLS_ENGEL_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 engel_expansion { | |
20 | public: | |
21 | engel_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 an Engel expansion."); | |
30 | } | |
31 | ||
32 | if(x==0) | |
33 | { | |
34 | throw std::domain_error("Zero does not have an Engel expansion."); | |
35 | } | |
36 | a_.reserve(64); | |
37 | // Let the error bound grow by 1 ULP/iteration. | |
38 | // I haven't done the error analysis to show that this is an expected rate of error growth, | |
39 | // but if you don't do this, you can easily get into an infinite loop. | |
40 | Real i = 1; | |
41 | Real computed = 0; | |
42 | Real term = 1; | |
43 | Real scale = std::numeric_limits<Real>::epsilon()*abs(x_)/2; | |
44 | Real u = x; | |
45 | while (abs(x_ - computed) > (i++)*scale) | |
46 | { | |
47 | Real recip = 1/u; | |
48 | Real ak = ceil(recip); | |
49 | a_.push_back(static_cast<Z>(ak)); | |
50 | u = u*ak - 1; | |
51 | if (u==0) | |
52 | { | |
53 | break; | |
54 | } | |
55 | term /= ak; | |
56 | computed += term; | |
57 | } | |
58 | ||
59 | for (size_t i = 1; i < a_.size(); ++i) | |
60 | { | |
61 | // Sanity check: This should only happen when wraparound occurs: | |
62 | if (a_[i] < a_[i-1]) | |
63 | { | |
64 | throw std::domain_error("The digits of an Engel expansion must form a non-decreasing sequence; consider increasing the wide of the integer type."); | |
65 | } | |
66 | // Watch out for saturating behavior: | |
67 | if (a_[i] == std::numeric_limits<Z>::max()) | |
68 | { | |
69 | throw std::domain_error("The integer type Z does not have enough width to hold the terms of the Engel expansion; please widen the type."); | |
70 | } | |
71 | } | |
72 | a_.shrink_to_fit(); | |
73 | } | |
74 | ||
75 | ||
76 | const std::vector<Z>& digits() const | |
77 | { | |
78 | return a_; | |
79 | } | |
80 | ||
81 | template<typename T, typename Z2> | |
82 | friend std::ostream& operator<<(std::ostream& out, engel_expansion<T, Z2>& eng); | |
83 | ||
84 | private: | |
85 | Real x_; | |
86 | std::vector<Z> a_; | |
87 | }; | |
88 | ||
89 | ||
90 | template<typename Real, typename Z2> | |
91 | std::ostream& operator<<(std::ostream& out, engel_expansion<Real, Z2>& engel) | |
92 | { | |
93 | constexpr const int p = std::numeric_limits<Real>::max_digits10; | |
94 | if constexpr (p == 2147483647) | |
95 | { | |
96 | out << std::setprecision(engel.x_.backend().precision()); | |
97 | } | |
98 | else | |
99 | { | |
100 | out << std::setprecision(p); | |
101 | } | |
102 | ||
103 | out << "{"; | |
104 | for (size_t i = 0; i < engel.a_.size() - 1; ++i) | |
105 | { | |
106 | out << engel.a_[i] << ", "; | |
107 | } | |
108 | out << engel.a_.back(); | |
109 | out << "}"; | |
110 | return out; | |
111 | } | |
112 | ||
113 | ||
114 | } | |
115 | #endif |