gum::Rational< GUM_SCALAR > Class Template Reference

Class template used to approximate decimal numbers by rationals. More...

#include <agrum/base/core/math/rational.h>

Static Public Member Functions
Real approximation by rational
static void	farey (int64_t &numerator, int64_t &denominator, const GUM_SCALAR &number, const int64_t &den_max=1000000L, const GUM_SCALAR &zero=1e-6)
	Find the rational close enough to a given ( decimal ) number in [-1,1] and whose denominator is not higher than a given integer number.
static void	continuedFracFirst (int64_t &numerator, int64_t &denominator, const GUM_SCALAR &number, const double &zero=1e-6)
	Find the first best rational approximation.
static void	continuedFracBest (int64_t &numerator, int64_t &denominator, const GUM_SCALAR &number, const int64_t &den_max=1000000)
	Find the best rational approximation.

Detailed Description

template<GUM_Numeric GUM_SCALAR>
class gum::Rational< GUM_SCALAR >

Class template used to approximate decimal numbers by rationals.

Template Parameters

GUM_SCALAR

The floating type ( float, double, long double ... ) of the number.

Definition at line 79 of file rational.h.

Member Function Documentation

continuedFracBest()

template<GUM_Numeric GUM_SCALAR>

void gum::Rational< GUM_SCALAR >::continuedFracBest ( int64_t & numerator, int64_t & denominator, const GUM_SCALAR & number, const int64_t & den_max = 1000000 )

static

Find the best rational approximation.

Not the first, to a given ( decimal) number ( ANY number ) and whose denominator is not higher than a given integer number.

In this case, we look for semi-convergents at the right of the last admissible convergent, if any. They are better approximations, but have higher denominators.

Parameters

numerator	The numerator of the rational.
denominator	The denominator of the rational.
number	The constant number we want to approximate using rationals.
den_max	The constant highest authorized denominator. 1000000 by default.

Definition at line 218 of file rational_tpl.h.

                                                                            {
    const GUM_SCALAR pnumber = (number > 0) ? number : -number;
 
    const uint64_t denMax = (uint64_t)den_max;   

    GUM_SCALAR rnumber = pnumber;

    std::vector< uint64_t > p({0, 1});
    std::vector< uint64_t > q({1, 0});

    std::vector< uint64_t > a;
 
    uint64_t p_tmp, q_tmp;
 
    uint64_t n;
    double   delta, delta_tmp;

    while (true) {
      a.push_back(std::lrint(std::floor(rnumber)));
 
      p_tmp = a.back() * p.back() + p[p.size() - 2];
      q_tmp = a.back() * q.back() + q[q.size() - 2];
 
      if (q_tmp > denMax || p_tmp > denMax) break;
 
      p.push_back(p_tmp);
      q.push_back(q_tmp);
 
      if (std::fabs(rnumber - a.back()) < 1e-6) break;
 
      rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
    }   
 
    if (a.size() < 2 || q.back() == denMax || p.back() == denMax) {
      numerator = p.back();
      if (number < 0) numerator = -numerator;
      denominator = q.back();
      return;
    }

    Idx i = Idx(p.size() - 1);

    for (n = a[i - 1] - 1; n >= (a[i - 1] + 2) / 2; --n) {
      p_tmp = n * p[i] + p[i - 1];
      q_tmp = n * q[i] + q[i - 1];
 
      if (q_tmp > denMax || p_tmp > denMax) continue;
 
      numerator = (int64_t)p_tmp;
      if (number < 0) numerator = -numerator;
      denominator = q_tmp;
      return;
    }   // end of for
 
    // Test n = a[i-1]/2
    n     = a[i - 1] / 2;
    p_tmp = n * p[i] + p[i - 1];
    q_tmp = n * q[i] + q[i - 1];
 
    delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
    delta     = std::fabs(pnumber - ((double)p[i]) / q[i]);
 
    if (delta_tmp < delta && q_tmp <= denMax && p_tmp <= denMax) {
      numerator = (int64_t)p_tmp;
      if (number < 0) numerator = -numerator;
      denominator = q_tmp;
    } else {
      numerator = (int64_t)p[i];
      if (number < 0) numerator = -numerator;
 
      denominator = q[i];
    }

  }

continuedFracFirst()

template<GUM_Numeric GUM_SCALAR>

void gum::Rational< GUM_SCALAR >::continuedFracFirst ( int64_t & numerator, int64_t & denominator, const GUM_SCALAR & number, const double & zero = 1e-6 )

static

Find the first best rational approximation.

end of farey func

The one with the smallest denominator such that no other rational with smaller denominator is a better approx, within precision zero to a given ( decimal ) number ( ANY number).

It gives the same answer than farey assuming zero is the same and den_max is infinite. Use this functions because you are sure to get an approx within zero of number.

We look at the semi-convergents left of the last admissible convergent, if any. They may be within the same precision and have a smaller denominator.

Parameters

numerator	The numerator of the rational.
denominator	The denominator of the rational.
number	The constant number we want to approximate using rationals.
zero	The positive value below which a number is considered zero. 1e-6 by default.

Definition at line 121 of file rational_tpl.h.

                                                                          {
    const GUM_SCALAR pnumber = (number > 0) ? number : -number;

    GUM_SCALAR rnumber = pnumber;

    std::vector< uint64_t > p({0, 1});
    std::vector< uint64_t > q({1, 0});

    std::vector< uint64_t > a;
 
    uint64_t p_tmp, q_tmp;
 
    uint64_t n;
    double   delta, delta_tmp;

    while (true) {
      a.push_back(std::lrint(std::floor(rnumber)));
      p.push_back(a.back() * p.back() + p[p.size() - 2]);
      q.push_back(a.back() * q.back() + q[q.size() - 2]);
 
      delta = std::fabs(pnumber - (GUM_SCALAR)p.back() / q.back());
 
      if (delta < zero) {
        numerator = (int64_t)p.back();
        if (number < 0) numerator = -numerator;
        denominator = q.back();
        break;
      }
 
      if (std::abs(rnumber - a.back()) < 1e-6) break;
 
      rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
    }   
 
    if (a.size() < 2) return;

    Idx i = Idx(p.size() - 2);
    // Test n = a[i-1]/2 ( when a[i-1] is even )
    n     = a[i - 1] / 2;
    p_tmp = n * p[i] + p[i - 1];
    q_tmp = n * q[i] + q[i - 1];
 
    delta     = std::fabs(pnumber - ((double)p[i]) / q[i]);
    delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
 
    if (delta < zero) {
      numerator = (int64_t)p[i];
      if (number < 0) numerator = -numerator;
      denominator = q[i];
      return;
    }
 
    if (delta_tmp < zero) {
      numerator = (int64_t)p_tmp;
      if (number < 0) numerator = -numerator;
      denominator = q_tmp;
      return;
    }
 
    // next semi-convergents until next convergent from smaller denominator to
    // bigger
    // denominator
    for (n = (a[i - 1] + 2) / 2; n < a[i - 1]; ++n) {
      p_tmp = n * p[i] + p[i - 1];
      q_tmp = n * q[i] + q[i - 1];
 
      delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
 
      if (delta_tmp < zero) {
        numerator = (int64_t)p_tmp;
        if (number < 0) numerator = -numerator;
        denominator = q_tmp;
        return;
      }
    }   

  }

Referenced by gum::credal::LRSWrapper< GUM_SCALAR >::_fill_().

Here is the caller graph for this function:

farey()

template<GUM_Numeric GUM_SCALAR>

void gum::Rational< GUM_SCALAR >::farey ( int64_t & numerator, int64_t & denominator, const GUM_SCALAR & number, const int64_t & den_max = 1000000L, const GUM_SCALAR & zero = 1e-6 )

static

Find the rational close enough to a given ( decimal ) number in [-1,1] and whose denominator is not higher than a given integer number.

Because of the double constraint on precision and size of the denominator, there is no guarantee on the precision of the approximation if den_max is low and zero is high. Prefer the use of continued fractions.

Parameters

numerator	The numerator of the rational.
denominator	The denominator of the rational.
number	The constant number we want to approximate using rationals.
den_max	The constant highest authorized denominator. 1000000 by default.
zero	The positive value below which a number is considered zero. 1e-6 by default.

Definition at line 58 of file rational_tpl.h.

                                                             {
    bool       isNegative = (number < 0) ? true : false;
    GUM_SCALAR pnumber    = (isNegative) ? -number : number;
 
    if (std::abs(pnumber - GUM_SCALAR(1.)) < zero) {
      numerator   = (isNegative) ? -1 : 1;
      denominator = 1;
      return;
    } else if (pnumber < zero) {
      numerator   = 0;
      denominator = 1;
      return;
    }
 
    int64_t a(0), b(1), c(1), d(1);
    double  mediant(0.0F);
 
    while (b <= den_max && d <= den_max) {
      // guard against int64_t overflow before computing a+c and b+d
      if (c > std::numeric_limits< int64_t >::max() - a
          || d > std::numeric_limits< int64_t >::max() - b)
        break;
      mediant = (GUM_SCALAR)(a + c) / (GUM_SCALAR)(b + d);
 
      if (std::fabs(pnumber - mediant) < zero) {
        if (b + d <= den_max) {
          numerator   = (isNegative) ? -(a + c) : (a + c);
          denominator = b + d;
          return;
        } else if (d > b) {
          numerator   = (isNegative) ? -c : c;
          denominator = d;
          return;
        } else {
          numerator   = (isNegative) ? -a : a;
          denominator = b;
          return;
        }
      } else if (pnumber > mediant) {
        a = a + c;
        b = b + d;
      } else {
        c = a + c;
        d = b + d;
      }
    }
 
    if (b > den_max) {
      numerator   = (isNegative) ? -c : c;
      denominator = d;
      return;
    } else {
      numerator   = (isNegative) ? -a : a;
      denominator = b;
      return;
    }
  }   

Referenced by gum::credal::CredalNet< GUM_SCALAR >::_H2Vlrs_().

Here is the caller graph for this function:

---

The documentation for this class was generated from the following files:

• agrum/base/core/math/rational.h
• agrum/base/core/math/rational_tpl.h