aGrUM 3.1.1
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-2026 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-2026 *
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
41#pragma once
42
43
49
50// To help IDE Parsers
51#include <agrum/agrum.h>
52
53#include <agrum/base/core/math/rational.h>
54
55namespace gum {
56
57 template < GUM_Numeric GUM_SCALAR >
58 void Rational< GUM_SCALAR >::farey(int64_t& numerator,
59 int64_t& denominator,
60 const GUM_SCALAR& number,
61 const int64_t& den_max,
62 const GUM_SCALAR& zero) {
63 bool isNegative = (number < 0) ? true : false;
64 GUM_SCALAR pnumber = (isNegative) ? -number : number;
65
66 if (std::abs(pnumber - GUM_SCALAR(1.)) < zero) {
67 numerator = (isNegative) ? -1 : 1;
68 denominator = 1;
69 return;
70 } else if (pnumber < zero) {
71 numerator = 0;
72 denominator = 1;
73 return;
74 }
75
76 int64_t a(0), b(1), c(1), d(1);
77 double mediant(0.0F);
78
79 while (b <= den_max && d <= den_max) {
80 // guard against int64_t overflow before computing a+c and b+d
81 if (c > std::numeric_limits< int64_t >::max() - a
82 || d > std::numeric_limits< int64_t >::max() - b)
83 break;
84 mediant = (GUM_SCALAR)(a + c) / (GUM_SCALAR)(b + d);
85
86 if (std::fabs(pnumber - mediant) < zero) {
87 if (b + d <= den_max) {
88 numerator = (isNegative) ? -(a + c) : (a + c);
89 denominator = b + d;
90 return;
91 } else if (d > b) {
92 numerator = (isNegative) ? -c : c;
93 denominator = d;
94 return;
95 } else {
96 numerator = (isNegative) ? -a : a;
97 denominator = b;
98 return;
99 }
100 } else if (pnumber > mediant) {
101 a = a + c;
102 b = b + d;
103 } else {
104 c = a + c;
105 d = b + d;
106 }
107 }
108
109 if (b > den_max) {
110 numerator = (isNegative) ? -c : c;
111 denominator = d;
112 return;
113 } else {
114 numerator = (isNegative) ? -a : a;
115 denominator = b;
116 return;
117 }
118 }
119
120 template < GUM_Numeric GUM_SCALAR >
121 void Rational< GUM_SCALAR >::continuedFracFirst(int64_t& numerator,
122 int64_t& denominator,
123 const GUM_SCALAR& number,
124 const double& zero) {
125 const GUM_SCALAR pnumber = (number > 0) ? number : -number;
126
128 GUM_SCALAR rnumber = pnumber;
129
131 std::vector< uint64_t > p({0, 1});
132 std::vector< uint64_t > q({1, 0});
133
135 std::vector< uint64_t > a;
136
137 uint64_t p_tmp, q_tmp;
138
139 uint64_t n;
140 double delta, delta_tmp;
141
147 while (true) {
148 a.push_back(std::lrint(std::floor(rnumber)));
149 p.push_back(a.back() * p.back() + p[p.size() - 2]);
150 q.push_back(a.back() * q.back() + q[q.size() - 2]);
151
152 delta = std::fabs(pnumber - (GUM_SCALAR)p.back() / q.back());
153
154 if (delta < zero) {
155 numerator = (int64_t)p.back();
156 if (number < 0) numerator = -numerator;
157 denominator = q.back();
158 break;
159 }
160
161 if (std::abs(rnumber - a.back()) < 1e-6) break;
162
163 rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
164 }
165
166 if (a.size() < 2) return;
167
171 Idx i = Idx(p.size() - 2);
175 // Test n = a[i-1]/2 ( when a[i-1] is even )
176 n = a[i - 1] / 2;
177 p_tmp = n * p[i] + p[i - 1];
178 q_tmp = n * q[i] + q[i - 1];
179
180 delta = std::fabs(pnumber - ((double)p[i]) / q[i]);
181 delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
182
183 if (delta < zero) {
184 numerator = (int64_t)p[i];
185 if (number < 0) numerator = -numerator;
186 denominator = q[i];
187 return;
188 }
189
190 if (delta_tmp < zero) {
191 numerator = (int64_t)p_tmp;
192 if (number < 0) numerator = -numerator;
193 denominator = q_tmp;
194 return;
195 }
196
197 // next semi-convergents until next convergent from smaller denominator to
198 // bigger
199 // denominator
200 for (n = (a[i - 1] + 2) / 2; n < a[i - 1]; ++n) {
201 p_tmp = n * p[i] + p[i - 1];
202 q_tmp = n * q[i] + q[i - 1];
203
204 delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
205
206 if (delta_tmp < zero) {
207 numerator = (int64_t)p_tmp;
208 if (number < 0) numerator = -numerator;
209 denominator = q_tmp;
210 return;
211 }
212 }
213
215 }
216
217 template < GUM_Numeric GUM_SCALAR >
218 void Rational< GUM_SCALAR >::continuedFracBest(int64_t& numerator,
219 int64_t& denominator,
220 const GUM_SCALAR& number,
221 const int64_t& den_max) {
222 const GUM_SCALAR pnumber = (number > 0) ? number : -number;
223
224 const uint64_t denMax = (uint64_t)den_max;
225
227 GUM_SCALAR rnumber = pnumber;
228
230 std::vector< uint64_t > p({0, 1});
231 std::vector< uint64_t > q({1, 0});
232
234 std::vector< uint64_t > a;
235
236 uint64_t p_tmp, q_tmp;
237
238 uint64_t n;
239 double delta, delta_tmp;
240
242 while (true) {
243 a.push_back(std::lrint(std::floor(rnumber)));
244
245 p_tmp = a.back() * p.back() + p[p.size() - 2];
246 q_tmp = a.back() * q.back() + q[q.size() - 2];
247
248 if (q_tmp > denMax || p_tmp > denMax) break;
249
250 p.push_back(p_tmp);
251 q.push_back(q_tmp);
252
253 if (std::fabs(rnumber - a.back()) < 1e-6) break;
254
255 rnumber = GUM_SCALAR(1.) / (rnumber - a.back());
256 }
257
258 if (a.size() < 2 || q.back() == denMax || p.back() == denMax) {
259 numerator = p.back();
260 if (number < 0) numerator = -numerator;
261 denominator = q.back();
262 return;
263 }
264
268 Idx i = Idx(p.size() - 1);
269
273 for (n = a[i - 1] - 1; n >= (a[i - 1] + 2) / 2; --n) {
274 p_tmp = n * p[i] + p[i - 1];
275 q_tmp = n * q[i] + q[i - 1];
276
277 if (q_tmp > denMax || p_tmp > denMax) continue;
278
279 numerator = (int64_t)p_tmp;
280 if (number < 0) numerator = -numerator;
281 denominator = q_tmp;
282 return;
283 } // end of for
284
285 // Test n = a[i-1]/2
286 n = a[i - 1] / 2;
287 p_tmp = n * p[i] + p[i - 1];
288 q_tmp = n * q[i] + q[i - 1];
289
290 delta_tmp = std::fabs(pnumber - ((double)p_tmp) / q_tmp);
291 delta = std::fabs(pnumber - ((double)p[i]) / q[i]);
292
293 if (delta_tmp < delta && q_tmp <= denMax && p_tmp <= denMax) {
294 numerator = (int64_t)p_tmp;
295 if (number < 0) numerator = -numerator;
296 denominator = q_tmp;
297 } else {
298 numerator = (int64_t)p[i];
299 if (number < 0) numerator = -numerator;
300
301 denominator = q[i];
302 }
303
305 }
306
307} // 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:218
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:121
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:58
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 3.1.1

© PHW&CG&others - 2022
DoXyGeN 1.16.1