aGrUM  0.13.2
rational_tpl.h
Go to the documentation of this file.
1 /***************************************************************************
2  * Copyright (C) 2005 by Pierre-Henri WUILLEMIN and Christophe GONZALES *
3  * {prenom.nom}_at_lip6.fr *
4  * *
5  * This program is free software; you can redistribute it and/or modify *
6  * it under the terms of the GNU General Public License as published by *
7  * the Free Software Foundation; either version 2 of the License, or *
8  * (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public License *
16  * along with this program; if not, write to the *
17  * Free Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ***************************************************************************/
20 
27 // To help IDE Parsers
28 #include <agrum/agrum.h>
30 
31 namespace gum {
32 
33  template < typename GUM_SCALAR >
34  void Rational< GUM_SCALAR >::farey(int64_t& numerator,
35  int64_t& denominator,
36  const GUM_SCALAR& number,
37  const int64_t& den_max,
38  const GUM_SCALAR& zero) {
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  }
91 
92  template < typename GUM_SCALAR >
94  int64_t& denominator,
95  const GUM_SCALAR& number,
96  const double& zero) {
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  }
188 
189  template < typename GUM_SCALAR >
191  int64_t& denominator,
192  const GUM_SCALAR& number,
193  const int64_t& den_max) {
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  }
279 
280 } // namespace gum
gum is the global namespace for all aGrUM entities
Definition: agrum.h:25
static void continuedFracBest(int64_t &numerator, int64_t &denominator, const GUM_SCALAR &number, const int64_t &den_max=1000000)
Find the best rational approximation.
Definition: rational_tpl.h:190
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 h...
Definition: rational_tpl.h:34
static void continuedFracFirst(int64_t &numerator, int64_t &denominator, const GUM_SCALAR &number, const double &zero=1e-6)
Find the first best rational approximation.
Definition: rational_tpl.h:93
Class template used to approximate decimal numbers by rationals.
unsigned long Idx
Type for indexes.
Definition: types.h:43