aGrUM  0.14.2
gum::Rational< GUM_SCALAR > Class Template Reference

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

#include <agrum/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 57 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 190 of file rational_tpl.h.

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

◆ 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 93 of file rational_tpl.h.

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

96  {
97  const GUM_SCALAR pnumber = (number > 0) ? number : -number;
98 
100  GUM_SCALAR rnumber = pnumber;
101 
103  std::vector< uint64_t > p({0, 1});
104  std::vector< uint64_t > q({1, 0});
105 
107  std::vector< uint64_t > a;
108 
109  uint64_t p_tmp, q_tmp;
110 
111  uint64_t n;
112  double delta, delta_tmp;
113 
119  while (true) {
120  a.push_back(std::lrint(std::floor(rnumber)));
121  p.push_back(a.back() * p.back() + p[p.size() - 2]);
122  q.push_back(a.back() * q.back() + q[q.size() - 2]);
123 
124  delta = std::fabs(pnumber - (GUM_SCALAR)p.back() / q.back());
125 
126  if (delta < zero) {
127  numerator = (int64_t)p.back();
128  if (number < 0) numerator = -numerator;
129  denominator = q.back();
130  break;
131  }
132 
133  if (std::abs(rnumber - a.back()) < 1e-6) break;
134 
135  rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
136  }
137 
138  if (a.size() < 2) return;
139 
143  Idx i = Idx(p.size() - 2);
147  // Test n = a[i-1]/2 ( when a[i-1] is even )
148  n = a[i - 1] / 2;
149  p_tmp = n * p[i] + p[i - 1];
150  q_tmp = n * q[i] + q[i - 1];
151 
152  delta = std::fabs(pnumber - ((double)p[i]) / q[i]);
153  delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
154 
155  if (delta < zero) {
156  numerator = (int64_t)p[i];
157  if (number < 0) numerator = -numerator;
158  denominator = q[i];
159  return;
160  }
161 
162  if (delta_tmp < zero) {
163  numerator = (int64_t)p_tmp;
164  if (number < 0) numerator = -numerator;
165  denominator = q_tmp;
166  return;
167  }
168 
169  // next semi-convergents until next convergent from smaller denominator to
170  // bigger
171  // denominator
172  for (n = (a[i - 1] + 2) / 2; n < a[i - 1]; ++n) {
173  p_tmp = n * p[i] + p[i - 1];
174  q_tmp = n * q[i] + q[i - 1];
175 
176  delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
177 
178  if (delta_tmp < zero) {
179  numerator = (int64_t)p_tmp;
180  if (number < 0) numerator = -numerator;
181  denominator = q_tmp;
182  return;
183  }
184  }
185 
187  }
Size Idx
Type for indexes.
Definition: types.h:50
+ 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 34 of file rational_tpl.h.

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

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

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