Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <cmath>
8 : #include <cstdint>
9 : #include <limits>
10 : #include <type_traits>
11 :
12 : #include "Utilities/ConstantExpressions.hpp"
13 : #include "Utilities/Gsl.hpp"
14 :
15 : /// Type trait to check if a type looks like a fraction (specifically,
16 : /// if it has numerator and denominator methods)
17 : /// @{
18 : template <typename T, typename = std::void_t<>>
19 1 : struct is_fraction : std::false_type {};
20 : template <typename T>
21 : struct is_fraction<T, std::void_t<decltype(std::declval<T>().numerator()),
22 : decltype(std::declval<T>().denominator())>>
23 : : std::true_type {};
24 :
25 : template <typename T>
26 0 : constexpr bool is_fraction_v = is_fraction<T>::value;
27 : /// @}
28 :
29 : /// \ingroup UtilitiesGroup
30 : /// \brief Compute the continued fraction representation of a number
31 : ///
32 : /// The template argument may be a fraction type or a floating point
33 : /// type. If the expansion is computed exactly then it will be the
34 : /// shorter of the expansions for the given value.
35 : template <class T>
36 1 : class ContinuedFraction {
37 : template <typename U, typename = typename is_fraction<U>::type>
38 0 : struct value_type_helper;
39 :
40 : template <typename U>
41 : struct value_type_helper<U, std::false_type> {
42 : using type = int64_t;
43 : };
44 :
45 : template <typename U>
46 : struct value_type_helper<U, std::true_type> {
47 : using type = std::decay_t<decltype(std::declval<U>().numerator())>;
48 : };
49 :
50 : public:
51 : /// The decayed return type of T::numerator() for a fraction,
52 : /// int64_t for other types.
53 1 : using value_type = typename value_type_helper<T>::type;
54 :
55 0 : explicit ContinuedFraction(const T& value)
56 : : term_(ifloor(value)),
57 : remainder_(value - term_),
58 : error_([this, &value]() {
59 : using std::abs;
60 : using std::max;
61 : // For any non-fundamental type, epsilon() returns a
62 : // default-constructed value, which should be zero for
63 : // fractions.
64 : return max(abs(value), abs(remainder_)) *
65 : std::numeric_limits<T>::epsilon();
66 : }()) {}
67 :
68 : /// Obtain the current element in the expansion
69 1 : value_type operator*() const { return term_; }
70 :
71 : /// Check if the expansion is incomplete. If T is a fraction type
72 : /// this will be true until the full exact representation of the
73 : /// fraction is produced. If T is a floating point type this will
74 : /// be true until the estimated numerical error is larger than the
75 : /// remaining error in the representation.
76 1 : explicit operator bool() const { return not done_; }
77 :
78 : /// Advance to the next element in the expansion
79 1 : ContinuedFraction& operator++() {
80 : // Terminate when remainder_ is consistent with zero.
81 : if (remainder_ == 0 or error_ > remainder_) {
82 : done_ = true;
83 : return *this;
84 : }
85 : remainder_ = 1 / remainder_;
86 : error_ *= square(remainder_);
87 : term_ = ifloor(remainder_);
88 : remainder_ -= term_;
89 : return *this;
90 : }
91 :
92 : private:
93 : template <typename U>
94 0 : static value_type ifloor(const U& x) {
95 : if constexpr (is_fraction_v<U>) {
96 : return static_cast<value_type>(x.numerator() / x.denominator());
97 : } else {
98 : return static_cast<value_type>(std::floor(x));
99 : }
100 : }
101 :
102 0 : value_type term_;
103 0 : T remainder_;
104 : // Estimate of error in the term.
105 0 : T error_;
106 0 : bool done_{false};
107 : };
108 :
109 : /// \ingroup UtilitiesGroup
110 : /// \brief Sum a continued fraction
111 : ///
112 : /// \tparam Fraction the result type, which must be a fraction type
113 : template <class Fraction>
114 1 : class ContinuedFractionSummer {
115 : public:
116 0 : using Term_t = std::decay_t<decltype(std::declval<Fraction>().numerator())>;
117 :
118 : /// Sum of the supplied continued fraction terms. If the exact sum
119 : /// cannot be represented by the `Fraction` type, a semiconvergent
120 : /// may be returned instead.
121 1 : Fraction value() const { return Fraction(numerator_, denominator_); }
122 :
123 : /// Insert a new term. Terms should be supplied from most to least
124 : /// significant.
125 1 : void insert(Term_t term) {
126 : if (overflowed_) {
127 : return;
128 : }
129 :
130 : Term_t new_numerator{};
131 : Term_t new_denominator{};
132 : if (__builtin_mul_overflow(term, numerator_, &new_numerator) or
133 : __builtin_add_overflow(new_numerator, prev_numerator_,
134 : &new_numerator) or
135 : __builtin_mul_overflow(term, denominator_, &new_denominator) or
136 : __builtin_add_overflow(new_denominator, prev_denominator_,
137 : &new_denominator)) {
138 : overflowed_ = true;
139 :
140 : // We can't add this term exactly, but we can try to return a
141 : // semiconvergent.
142 : using std::abs;
143 : const Term_t representable_semiconvergent_index =
144 : std::min((std::numeric_limits<Term_t>::max() - abs(prev_numerator_)) /
145 : abs(numerator_),
146 : (std::numeric_limits<Term_t>::max() - prev_denominator_) /
147 : denominator_);
148 : // There is an edge case where the two sides of this inequality
149 : // are equal, where we don't have enough information to
150 : // determine whether the semiconvergent approximation is better
151 : // than the current convergent. We assume it's not.
152 : if (2 * representable_semiconvergent_index <= term) {
153 : return;
154 : }
155 :
156 : new_numerator =
157 : representable_semiconvergent_index * numerator_ + prev_numerator_;
158 : new_denominator =
159 : representable_semiconvergent_index * denominator_ + prev_denominator_;
160 : }
161 :
162 : prev_numerator_ = numerator_;
163 : numerator_ = new_numerator;
164 : prev_denominator_ = denominator_;
165 : denominator_ = new_denominator;
166 : }
167 :
168 : private:
169 0 : Term_t numerator_{1};
170 0 : Term_t denominator_{0};
171 0 : Term_t prev_numerator_{0};
172 0 : Term_t prev_denominator_{1};
173 0 : bool overflowed_{false};
174 : };
175 :
176 : /// \ingroup UtilitiesGroup
177 : /// \brief Find the fraction in the supplied interval with the
178 : /// smallest denominator
179 : ///
180 : /// The endpoints are considered to be in the interval. The order of
181 : /// the arguments is not significant. The answer is unique as long as
182 : /// the interval has length less than 1; for longer intervals, an
183 : /// integer in the range will be returned.
184 : template <typename Fraction, typename T1, typename T2>
185 1 : Fraction simplest_fraction_in_interval(const T1& end1, const T2& end2) {
186 : ContinuedFractionSummer<Fraction> result;
187 : using Term_t = typename decltype(result)::Term_t;
188 : ContinuedFraction<T1> cf1(end1);
189 : ContinuedFraction<T2> cf2(end2);
190 : using InputTerm_t = std::common_type_t<typename decltype(cf1)::value_type,
191 : typename decltype(cf2)::value_type>;
192 : InputTerm_t term1 = *cf1;
193 : InputTerm_t term2 = *cf2;
194 : for (; cf1 and cf2; term1 = *++cf1, term2 = *++cf2) {
195 : if (term1 != term2) {
196 : if (++cf1) {
197 : ++term1;
198 : }
199 : if (++cf2) {
200 : ++term2;
201 : }
202 : using std::min;
203 : result.insert(gsl::narrow<Term_t>(min(term1, term2)));
204 : return result.value();
205 : }
206 : result.insert(gsl::narrow<Term_t>(term1));
207 : }
208 : return result.value();
209 : }
|
