SpECTRE Documentation Coverage Report
Current view: top level - Utilities - FractionUtilities.hpp Hit Total Coverage
Commit: b5f497991094937944b0a3f519166bb54739d08a Lines: 10 25 40.0 %
Date: 2024-03-28 18:20:13
Legend: Lines: hit not hit

          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
     119           1 :   Fraction value() const { return Fraction(numerator_, denominator_); }
     120             : 
     121             :   /// Insert a new term.  Terms should be supplied from most to least
     122             :   /// significant.
     123           1 :   void insert(Term_t term) {
     124             :     if (overflowed_) {
     125             :       return;
     126             :     }
     127             : 
     128             :     Term_t new_numerator{};
     129             :     if (__builtin_mul_overflow(term, numerator_, &new_numerator) or
     130             :         __builtin_add_overflow(new_numerator, prev_numerator_,
     131             :                                &new_numerator)) {
     132             :       overflowed_ = true;
     133             :       return;
     134             :     }
     135             : 
     136             :     Term_t new_denominator{};
     137             :     if (__builtin_mul_overflow(term, denominator_, &new_denominator) or
     138             :         __builtin_add_overflow(new_denominator, prev_denominator_,
     139             :                                &new_denominator)) {
     140             :       overflowed_ = true;
     141             :       return;
     142             :     }
     143             : 
     144             :     prev_numerator_ = numerator_;
     145             :     numerator_ = new_numerator;
     146             :     prev_denominator_ = denominator_;
     147             :     denominator_ = new_denominator;
     148             :   }
     149             : 
     150             :  private:
     151           0 :   Term_t numerator_{1};
     152           0 :   Term_t denominator_{0};
     153           0 :   Term_t prev_numerator_{0};
     154           0 :   Term_t prev_denominator_{1};
     155           0 :   bool overflowed_{false};
     156             : };
     157             : 
     158             : /// \ingroup UtilitiesGroup
     159             : /// \brief Find the fraction in the supplied interval with the
     160             : /// smallest denominator
     161             : ///
     162             : /// The endpoints are considered to be in the interval.  The order of
     163             : /// the arguments is not significant.  The answer is unique as long as
     164             : /// the interval has length less than 1; for longer intervals, an
     165             : /// integer in the range will be returned.
     166             : template <typename Fraction, typename T1, typename T2>
     167           1 : Fraction simplest_fraction_in_interval(const T1& end1, const T2& end2) {
     168             :   ContinuedFractionSummer<Fraction> result;
     169             :   using Term_t = typename decltype(result)::Term_t;
     170             :   ContinuedFraction<T1> cf1(end1);
     171             :   ContinuedFraction<T2> cf2(end2);
     172             :   using InputTerm_t = std::common_type_t<typename decltype(cf1)::value_type,
     173             :                                          typename decltype(cf2)::value_type>;
     174             :   InputTerm_t term1 = *cf1;
     175             :   InputTerm_t term2 = *cf2;
     176             :   for (; cf1 and cf2; term1 = *++cf1, term2 = *++cf2) {
     177             :     if (term1 != term2) {
     178             :       if (++cf1) {
     179             :         ++term1;
     180             :       }
     181             :       if (++cf2) {
     182             :         ++term2;
     183             :       }
     184             :       using std::min;
     185             :       result.insert(gsl::narrow<Term_t>(min(term1, term2)));
     186             :       return result.value();
     187             :     }
     188             :     result.insert(gsl::narrow<Term_t>(term1));
     189             :   }
     190             :   return result.value();
     191             : }

Generated by: LCOV version 1.14