aGrUM 2.3.2
a C++ library for (probabilistic) graphical models
rational_tpl.h
Go to the documentation of this file.
1/****************************************************************************
2 * This file is part of the aGrUM/pyAgrum library. *
3 * *
4 * Copyright (c) 2005-2025 by *
5 * - Pierre-Henri WUILLEMIN(_at_LIP6) *
6 * - Christophe GONZALES(_at_AMU) *
7 * *
8 * The aGrUM/pyAgrum library is free software; you can redistribute it *
9 * and/or modify it under the terms of either : *
10 * *
11 * - the GNU Lesser General Public License as published by *
12 * the Free Software Foundation, either version 3 of the License, *
13 * or (at your option) any later version, *
14 * - the MIT license (MIT), *
15 * - or both in dual license, as here. *
16 * *
17 * (see https://agrum.gitlab.io/articles/dual-licenses-lgplv3mit.html) *
18 * *
19 * This aGrUM/pyAgrum library is distributed in the hope that it will be *
20 * useful, but WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, *
21 * INCLUDING BUT NOT LIMITED TO THE WARRANTIES MERCHANTABILITY or FITNESS *
22 * FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE *
23 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
25 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR *
26 * OTHER DEALINGS IN THE SOFTWARE. *
27 * *
28 * See LICENCES for more details. *
29 * *
30 * SPDX-FileCopyrightText: Copyright 2005-2025 *
31 * - Pierre-Henri WUILLEMIN(_at_LIP6) *
32 * - Christophe GONZALES(_at_AMU) *
33 * SPDX-License-Identifier: LGPL-3.0-or-later OR MIT *
34 * *
35 * Contact : info_at_agrum_dot_org *
36 * homepage : http://agrum.gitlab.io *
37 * gitlab : https://gitlab.com/agrumery/agrum *
38 * *
39 ****************************************************************************/
40#pragma once
41
42
48
49// To help IDE Parsers
50#include <agrum/agrum.h>
51
52#include <agrum/base/core/math/rational.h>
53
54namespace gum {
55
56 template < typename GUM_SCALAR >
57 void Rational< GUM_SCALAR >::farey(int64_t& numerator,
58 int64_t& denominator,
59 const GUM_SCALAR& number,
60 const int64_t& den_max,
61 const GUM_SCALAR& zero) {
62 bool isNegative = (number < 0) ? true : false;
63 GUM_SCALAR pnumber = (isNegative) ? -number : number;
64
65 if (std::abs(pnumber - GUM_SCALAR(1.)) < zero) {
66 numerator = (isNegative) ? -1 : 1;
67 denominator = 1;
68 return;
69 } else if (pnumber < zero) {
70 numerator = 0;
71 denominator = 1;
72 return;
73 }
74
75 int64_t a(0), b(1), c(1), d(1);
76 double mediant(0.0F);
77
78 while (b <= den_max && d <= den_max) {
79 mediant = (GUM_SCALAR)(a + c) / (GUM_SCALAR)(b + d);
80
81 if (std::fabs(pnumber - mediant) < zero) {
82 if (b + d <= den_max) {
83 numerator = (isNegative) ? -(a + c) : (a + c);
84 denominator = b + d;
85 return;
86 } else if (d > b) {
87 numerator = (isNegative) ? -c : c;
88 denominator = d;
89 return;
90 } else {
91 numerator = (isNegative) ? -a : a;
92 denominator = b;
93 return;
94 }
95 } else if (pnumber > mediant) {
96 a = a + c;
97 b = b + d;
98 } else {
99 c = a + c;
100 d = b + d;
101 }
102 }
103
104 if (b > den_max) {
105 numerator = (isNegative) ? -c : c;
106 denominator = d;
107 return;
108 } else {
109 numerator = (isNegative) ? -a : a;
110 denominator = b;
111 return;
112 }
113 }
114
115 template < typename GUM_SCALAR >
116 void Rational< GUM_SCALAR >::continuedFracFirst(int64_t& numerator,
117 int64_t& denominator,
118 const GUM_SCALAR& number,
119 const double& zero) {
120 const GUM_SCALAR pnumber = (number > 0) ? number : -number;
121
123 GUM_SCALAR rnumber = pnumber;
124
126 std::vector< uint64_t > p({0, 1});
127 std::vector< uint64_t > q({1, 0});
128
130 std::vector< uint64_t > a;
131
132 uint64_t p_tmp, q_tmp;
133
134 uint64_t n;
135 double delta, delta_tmp;
136
142 while (true) {
143 a.push_back(std::lrint(std::floor(rnumber)));
144 p.push_back(a.back() * p.back() + p[p.size() - 2]);
145 q.push_back(a.back() * q.back() + q[q.size() - 2]);
146
147 delta = std::fabs(pnumber - (GUM_SCALAR)p.back() / q.back());
148
149 if (delta < zero) {
150 numerator = (int64_t)p.back();
151 if (number < 0) numerator = -numerator;
152 denominator = q.back();
153 break;
154 }
155
156 if (std::abs(rnumber - a.back()) < 1e-6) break;
157
158 rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
159 }
160
161 if (a.size() < 2) return;
162
166 Idx i = Idx(p.size() - 2);
170 // Test n = a[i-1]/2 ( when a[i-1] is even )
171 n = a[i - 1] / 2;
172 p_tmp = n * p[i] + p[i - 1];
173 q_tmp = n * q[i] + q[i - 1];
174
175 delta = std::fabs(pnumber - ((double)p[i]) / q[i]);
176 delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
177
178 if (delta < zero) {
179 numerator = (int64_t)p[i];
180 if (number < 0) numerator = -numerator;
181 denominator = q[i];
182 return;
183 }
184
185 if (delta_tmp < zero) {
186 numerator = (int64_t)p_tmp;
187 if (number < 0) numerator = -numerator;
188 denominator = q_tmp;
189 return;
190 }
191
192 // next semi-convergents until next convergent from smaller denominator to
193 // bigger
194 // denominator
195 for (n = (a[i - 1] + 2) / 2; n < a[i - 1]; ++n) {
196 p_tmp = n * p[i] + p[i - 1];
197 q_tmp = n * q[i] + q[i - 1];
198
199 delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
200
201 if (delta_tmp < zero) {
202 numerator = (int64_t)p_tmp;
203 if (number < 0) numerator = -numerator;
204 denominator = q_tmp;
205 return;
206 }
207 }
208
210 }
211
212 template < typename GUM_SCALAR >
213 void Rational< GUM_SCALAR >::continuedFracBest(int64_t& numerator,
214 int64_t& denominator,
215 const GUM_SCALAR& number,
216 const int64_t& den_max) {
217 const GUM_SCALAR pnumber = (number > 0) ? number : -number;
218
219 const uint64_t denMax = (uint64_t)den_max;
220
222 GUM_SCALAR rnumber = pnumber;
223
225 std::vector< uint64_t > p({0, 1});
226 std::vector< uint64_t > q({1, 0});
227
229 std::vector< uint64_t > a;
230
231 uint64_t p_tmp, q_tmp;
232
233 uint64_t n;
234 double delta, delta_tmp;
235
237 while (true) {
238 a.push_back(std::lrint(std::floor(rnumber)));
239
240 p_tmp = a.back() * p.back() + p[p.size() - 2];
241 q_tmp = a.back() * q.back() + q[q.size() - 2];
242
243 if (q_tmp > denMax || p_tmp > denMax) break;
244
245 p.push_back(p_tmp);
246 q.push_back(q_tmp);
247
248 if (std::fabs(rnumber - a.back()) < 1e-6) break;
249
250 rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
251 }
252
253 if (a.size() < 2 || q.back() == denMax || p.back() == denMax) {
254 numerator = p.back();
255 if (number < 0) numerator = -numerator;
256 denominator = q.back();
257 return;
258 }
259
263 Idx i = Idx(p.size() - 1);
264
268 for (n = a[i - 1] - 1; n >= (a[i - 1] + 2) / 2; --n) {
269 p_tmp = n * p[i] + p[i - 1];
270 q_tmp = n * q[i] + q[i - 1];
271
272 if (q_tmp > denMax || p_tmp > denMax) continue;
273
274 numerator = (int64_t)p_tmp;
275 if (number < 0) numerator = -numerator;
276 denominator = q_tmp;
277 return;
278 } // end of for
279
280 // Test n = a[i-1]/2
281 n = a[i - 1] / 2;
282 p_tmp = n * p[i] + p[i - 1];
283 q_tmp = n * q[i] + q[i - 1];
284
285 delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
286 delta = std::fabs(pnumber - ((double)p[i]) / q[i]);
287
288 if (delta_tmp < delta && q_tmp <= denMax && p_tmp <= denMax) {
289 numerator = (int64_t)p_tmp;
290 if (number < 0) numerator = -numerator;
291 denominator = q_tmp;
292 } else {
293 numerator = (int64_t)p[i];
294 if (number < 0) numerator = -numerator;
295
296 denominator = q[i];
297 }
298
300 }
301
302} // namespace gum
gum::Rational::continuedFracBest
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:213
gum::Rational::continuedFracFirst
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:116
gum::Rational::farey
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:57
gum::Idx
Size Idx
Type for indexes.
Definition types.h:79
gum
gum is the global namespace for all aGrUM entities
Definition agrum.h:46
rational.h
Class template used to approximate decimal numbers by rationals.
aGrUM 2.3.2

© PHW&CG&others - 2022
DoXyGeN 1.16.1