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             : }