aGrUM  0.18.1
a C++ library for (probabilistic) graphical models
gum::Rational< GUM_SCALAR > Class Template Reference

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

#include <agrum/tools/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. More...
 
static void continuedFracFirst (int64_t &numerator, int64_t &denominator, const GUM_SCALAR &number, const double &zero=1e-6)
 Find the first best rational approximation. More...
 
static void continuedFracBest (int64_t &numerator, int64_t &denominator, const GUM_SCALAR &number, const int64_t &den_max=1000000)
 Find the best rational approximation. More...
 

Detailed Description

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

Class template used to approximate decimal numbers by rationals.

Template Parameters
GUM_SCALARThe floating type ( float, double, long double ... ) of the number.

Definition at line 59 of file rational.h.

Member Function Documentation

◆ continuedFracBest()

template<typename 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
numeratorThe numerator of the rational.
denominatorThe denominator of the rational.
numberThe constant number we want to approximate using rationals.
den_maxThe constant highest authorized denominator. 1000000 by default.

Definition at line 192 of file rational_tpl.h.

195  {
196  const GUM_SCALAR pnumber = (number > 0) ? number : -number;
197 
198  const uint64_t denMax =
199  (uint64_t)den_max;
200 
202  GUM_SCALAR rnumber = pnumber;
203 
205  std::vector< uint64_t > p({0, 1});
206  std::vector< uint64_t > q({1, 0});
207 
209  std::vector< uint64_t > a;
210 
211  uint64_t p_tmp, q_tmp;
212 
213  uint64_t n;
214  double delta, delta_tmp;
215 
217  while (true) {
218  a.push_back(std::lrint(std::floor(rnumber)));
219 
220  p_tmp = a.back() * p.back() + p[p.size() - 2];
221  q_tmp = a.back() * q.back() + q[q.size() - 2];
222 
223  if (q_tmp > denMax || p_tmp > denMax) break;
224 
225  p.push_back(p_tmp);
226  q.push_back(q_tmp);
227 
228  if (std::fabs(rnumber - a.back()) < 1e-6) break;
229 
230  rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
231  }
232 
233  if (a.size() < 2 || q.back() == denMax || p.back() == denMax) {
234  numerator = p.back();
235  if (number < 0) numerator = -numerator;
236  denominator = q.back();
237  return;
238  }
239 
243  Idx i = Idx(p.size() - 1);
244 
248  for (n = a[i - 1] - 1; n >= (a[i - 1] + 2) / 2; --n) {
249  p_tmp = n * p[i] + p[i - 1];
250  q_tmp = n * q[i] + q[i - 1];
251 
252  if (q_tmp > denMax || p_tmp > denMax) continue;
253 
254  numerator = (int64_t)p_tmp;
255  if (number < 0) numerator = -numerator;
256  denominator = q_tmp;
257  return;
258  } // end of for
259 
260  // Test n = a[i-1]/2
261  n = a[i - 1] / 2;
262  p_tmp = n * p[i] + p[i - 1];
263  q_tmp = n * q[i] + q[i - 1];
264 
265  delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
266  delta = std::fabs(pnumber - ((double)p[i]) / q[i]);
267 
268  if (delta_tmp < delta && q_tmp <= denMax && p_tmp <= denMax) {
269  numerator = (int64_t)p_tmp;
270  if (number < 0) numerator = -numerator;
271  denominator = q_tmp;
272  } else {
273  numerator = (int64_t)p[i];
274  if (number < 0) numerator = -numerator;
275 
276  denominator = q[i];
277  }
278 
280  }
Size Idx
Type for indexes.
Definition: types.h:53

◆ continuedFracFirst()

template<typename 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
numeratorThe numerator of the rational.
denominatorThe denominator of the rational.
numberThe constant number we want to approximate using rationals.
zeroThe positive value below which a number is considered zero. 1e-6 by default.

Definition at line 95 of file rational_tpl.h.

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

98  {
99  const GUM_SCALAR pnumber = (number > 0) ? number : -number;
100 
102  GUM_SCALAR rnumber = pnumber;
103 
105  std::vector< uint64_t > p({0, 1});
106  std::vector< uint64_t > q({1, 0});
107 
109  std::vector< uint64_t > a;
110 
111  uint64_t p_tmp, q_tmp;
112 
113  uint64_t n;
114  double delta, delta_tmp;
115 
121  while (true) {
122  a.push_back(std::lrint(std::floor(rnumber)));
123  p.push_back(a.back() * p.back() + p[p.size() - 2]);
124  q.push_back(a.back() * q.back() + q[q.size() - 2]);
125 
126  delta = std::fabs(pnumber - (GUM_SCALAR)p.back() / q.back());
127 
128  if (delta < zero) {
129  numerator = (int64_t)p.back();
130  if (number < 0) numerator = -numerator;
131  denominator = q.back();
132  break;
133  }
134 
135  if (std::abs(rnumber - a.back()) < 1e-6) break;
136 
137  rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
138  }
139 
140  if (a.size() < 2) return;
141 
145  Idx i = Idx(p.size() - 2);
149  // Test n = a[i-1]/2 ( when a[i-1] is even )
150  n = a[i - 1] / 2;
151  p_tmp = n * p[i] + p[i - 1];
152  q_tmp = n * q[i] + q[i - 1];
153 
154  delta = std::fabs(pnumber - ((double)p[i]) / q[i]);
155  delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
156 
157  if (delta < zero) {
158  numerator = (int64_t)p[i];
159  if (number < 0) numerator = -numerator;
160  denominator = q[i];
161  return;
162  }
163 
164  if (delta_tmp < zero) {
165  numerator = (int64_t)p_tmp;
166  if (number < 0) numerator = -numerator;
167  denominator = q_tmp;
168  return;
169  }
170 
171  // next semi-convergents until next convergent from smaller denominator to
172  // bigger
173  // denominator
174  for (n = (a[i - 1] + 2) / 2; n < a[i - 1]; ++n) {
175  p_tmp = n * p[i] + p[i - 1];
176  q_tmp = n * q[i] + q[i - 1];
177 
178  delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
179 
180  if (delta_tmp < zero) {
181  numerator = (int64_t)p_tmp;
182  if (number < 0) numerator = -numerator;
183  denominator = q_tmp;
184  return;
185  }
186  }
187 
189  }
Size Idx
Type for indexes.
Definition: types.h:53
+ Here is the caller graph for this function:

◆ farey()

template<typename 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
numeratorThe numerator of the rational.
denominatorThe denominator of the rational.
numberThe constant number we want to approximate using rationals.
den_maxThe constant highest authorized denominator. 1000000 by default.
zeroThe positive value below which a number is considered zero. 1e-6 by default.

Definition at line 36 of file rational_tpl.h.

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

40  {
41  bool isNegative = (number < 0) ? true : false;
42  GUM_SCALAR pnumber = (isNegative) ? -number : number;
43 
44  if (std::abs(pnumber - GUM_SCALAR(1.)) < zero) {
45  numerator = (isNegative) ? -1 : 1;
46  denominator = 1;
47  return;
48  } else if (pnumber < zero) {
49  numerator = 0;
50  denominator = 1;
51  return;
52  }
53 
54  int64_t a(0), b(1), c(1), d(1);
55  double mediant(0.0F);
56 
57  while (b <= den_max && d <= den_max) {
58  mediant = (GUM_SCALAR)(a + c) / (GUM_SCALAR)(b + d);
59 
60  if (std::fabs(pnumber - mediant) < zero) {
61  if (b + d <= den_max) {
62  numerator = (isNegative) ? -(a + c) : (a + c);
63  denominator = b + d;
64  return;
65  } else if (d > b) {
66  numerator = (isNegative) ? -c : c;
67  denominator = d;
68  return;
69  } else {
70  numerator = (isNegative) ? -a : a;
71  denominator = b;
72  return;
73  }
74  } else if (pnumber > mediant) {
75  a = a + c;
76  b = b + d;
77  } else {
78  c = a + c;
79  d = b + d;
80  }
81  }
82 
83  if (b > den_max) {
84  numerator = (isNegative) ? -c : c;
85  denominator = d;
86  return;
87  } else {
88  numerator = (isNegative) ? -a : a;
89  denominator = b;
90  return;
91  }
92  }
+ Here is the caller graph for this function:

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