FractionUtilities.hpp
2 // See LICENSE.txt for details.
3
4 #pragma once
5
6 #include <cmath>
7 #include <cstdint>
8 #include <limits>
9 #include <type_traits>
10 #include <utility>
11
12 #include "ErrorHandling/Assert.hpp"
14 #include "Utilities/Gsl.hpp"
15 #include "Utilities/TypeTraits.hpp"
16
17 /// Type trait to check if a type looks like a fraction (specifically,
18 /// if it has numerator and denominator methods)
19 //@{
20 template <typename T, typename = cpp17::void_t<>>
22 template <typename T>
23 struct is_fraction<T, cpp17::void_t<decltype(std::declval<T>().numerator()),
24  decltype(std::declval<T>().denominator())>>
25  : std::true_type {};
26
27 template <typename T>
28 constexpr bool is_fraction_v = is_fraction<T>::value;
29 //@}
30
31 /// \ingroup UtilitiesGroup
32 /// \brief Compute the continued fraction representation of a number
33 ///
34 /// The template argument may be a fraction type or a floating point
35 /// type. If the expansion is computed exactly then it will be the
36 /// shorter of the expansions for the given value.
37 template <class T>
39  template <typename U, typename = typename is_fraction<U>::type>
40  struct value_type_helper;
41
42  template <typename U>
43  struct value_type_helper<U, std::false_type> {
44  using type = int64_t;
45  };
46
47  template <typename U>
48  struct value_type_helper<U, std::true_type> {
49  using type = std::decay_t<decltype(std::declval<U>().numerator())>;
50  };
51
52  public:
53  /// The decayed return type of T::numerator() for a fraction,
54  /// int64_t for other types.
55  using value_type = typename value_type_helper<T>::type;
56
57  explicit ContinuedFraction(const T& value) noexcept
58  : term_(ifloor(value)), remainder_(value - term_) {}
59
60  /// Obtain the current element in the expansion
61  value_type operator*() const noexcept { return term_; }
62
63  /// Check if the expansion is incomplete. If T is a fraction type
64  /// this will be true until the full exact representation of the
65  /// fraction is produced. If T is a floating point type this will
66  /// be true until the estimated numerical error is larger than the
67  /// remaining error in the representation.
68  explicit operator bool() const noexcept { return not done_; }
69
70  /// Advance to the next element in the expansion
72  if (remainder_ == 0 or (error_ /= square(remainder_)) > 1) {
73  done_ = true;
74  return *this;
75  }
76  remainder_ = 1 / remainder_;
77  term_ = ifloor(remainder_);
78  remainder_ -= term_;
79  return *this;
80  }
81
82  private:
83  template <typename U, Requires<is_fraction_v<U>> = nullptr>
84  static value_type ifloor(const U& x) noexcept {
85  return static_cast<value_type>(x.numerator() / x.denominator());
86  }
87
88  template <typename U, Requires<not is_fraction_v<U>> = nullptr>
89  static value_type ifloor(const U& x) noexcept {
90  return static_cast<value_type>(std::floor(x));
91  }
92
93  value_type term_;
94  T remainder_;
95  // Estimate of error in the term. For any non-fundamental type
96  // epsilon() returns a default-constructed value, which should be
97  // zero for fractions.
99  bool done_{false};
100 };
101
102 /// \ingroup UtilitiesGroup
103 /// \brief Sum a continued fraction
104 ///
105 /// \tparam Fraction the result type, which must be a fraction type
106 template <class Fraction>
108  public:
110
111  /// Sum of the supplied continued fraction terms
112  Fraction value() const noexcept { return Fraction(numerator_, denominator_); }
113
114  /// Insert a new term. Terms should be supplied from most to least
115  /// significant.
116  void insert(Term_t term) noexcept {
117  const auto new_numerator = term * numerator_ + prev_numerator_;
118  prev_numerator_ = numerator_;
119  numerator_ = new_numerator;
120
121  const auto new_denominator = term * denominator_ + prev_denominator_;
122  prev_denominator_ = denominator_;
123  denominator_ = new_denominator;
124  }
125
126  private:
127  Term_t numerator_{1};
128  Term_t denominator_{0};
129  Term_t prev_numerator_{0};
130  Term_t prev_denominator_{1};
131 };
132
133 /// \ingroup UtilitiesGroup
134 /// \brief Find the fraction in the supplied interval with the
135 /// smallest denominator
136 ///
137 /// The endpoints are considered to be in the interval. The order of
138 /// the arguments is not significant. The answer is unique as long as
139 /// the interval has length less than 1; for longer intervals, an
140 /// integer in the range will be returned.
141 template <typename Fraction, typename T1, typename T2>
142 Fraction simplest_fraction_in_interval(const T1& end1,
143  const T2& end2) noexcept {
145  using Term_t = typename decltype(result)::Term_t;
146  ContinuedFraction<T1> cf1(end1);
147  ContinuedFraction<T2> cf2(end2);
148  using InputTerm_t = std::common_type_t<typename decltype(cf1)::value_type,
149  typename decltype(cf2)::value_type>;
150  InputTerm_t term1 = *cf1;
151  InputTerm_t term2 = *cf2;
152  for (; cf1 and cf2; term1 = *++cf1, term2 = *++cf2) {
153  if (term1 != term2) {
154  if (++cf1) {
155  ++term1;
156  }
157  if (++cf2) {
158  ++term2;
159  }
160  result.insert(gsl::narrow<Term_t>(std::min(term1, term2)));
161  return result.value();
162  }
163  result.insert(gsl::narrow<Term_t>(term1));
164  }
165  return result.value();
166 }
ContinuedFraction & operator++() noexcept
Advance to the next element in the expansion.
Definition: FractionUtilities.hpp:71
void void_t
Given a set of types, returns void
Definition: TypeTraits.hpp:214
typename value_type_helper< T >::type value_type
The decayed return type of T::numerator() for a fraction, int64_t for other types.
Definition: FractionUtilities.hpp:55
T epsilon(T... args)
Type trait to check if a type looks like a fraction (specifically, if it has numerator and denominato...
Definition: FractionUtilities.hpp:21
Define simple functions for constant expressions.
void insert(Term_t term) noexcept
Insert a new term. Terms should be supplied from most to least significant.
Definition: FractionUtilities.hpp:116
Compute the continued fraction representation of a number.
Definition: FractionUtilities.hpp:38
Sum a continued fraction.
Definition: FractionUtilities.hpp:107
Defines macro ASSERT.
value_type operator*() const noexcept
Obtain the current element in the expansion.
Definition: FractionUtilities.hpp:61
Defines functions and classes from the GSL.
Fraction value() const noexcept
Sum of the supplied continued fraction terms.
Definition: FractionUtilities.hpp:112
Defines type traits, some of which are future STL type_traits header.
Fraction simplest_fraction_in_interval(const T1 &end1, const T2 &end2) noexcept
Find the fraction in the supplied interval with the smallest denominator.
Definition: FractionUtilities.hpp:142
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:52