 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_SCALAR The floating type ( float, double, long double ... ) of the number.

Definition at line 59 of file rational.h.

## ◆ 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
 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 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
 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 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
 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 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: