Branch data Line data Source code
1 : : /* SPDX-License-Identifier: BSD-3-Clause
2 : : * Copyright(c) 2010-2014 Intel Corporation
3 : : */
4 : :
5 : : #include <stdlib.h>
6 : :
7 : : #include <eal_export.h>
8 : : #include "rte_approx.h"
9 : :
10 : : /*
11 : : * Based on paper "Approximating Rational Numbers by Fractions" by Michal
12 : : * Forisek forisek@dcs.fmph.uniba.sk
13 : : *
14 : : * Given a rational number alpha with 0 < alpha < 1 and a precision d, the goal
15 : : * is to find positive integers p, q such that alpha - d < p/q < alpha + d, and
16 : : * q is minimal.
17 : : *
18 : : * http://people.ksp.sk/~misof/publications/2007approx.pdf
19 : : */
20 : :
21 : : /* fraction comparison: compare (a/b) and (c/d) */
22 : : static inline uint32_t
23 : : less(uint32_t a, uint32_t b, uint32_t c, uint32_t d)
24 : : {
25 : 0 : return a*d < b*c;
26 : : }
27 : :
28 : : static inline uint32_t
29 : : less_or_equal(uint32_t a, uint32_t b, uint32_t c, uint32_t d)
30 : : {
31 : 0 : return a*d <= b*c;
32 : : }
33 : :
34 : : /* check whether a/b is a valid approximation */
35 : : static inline uint32_t
36 : : matches(uint32_t a, uint32_t b,
37 : : uint32_t alpha_num, uint32_t d_num, uint32_t denum)
38 : : {
39 : 0 : if (less_or_equal(a, b, alpha_num - d_num, denum))
40 : : return 0;
41 : :
42 [ # # # # : 0 : if (less(a ,b, alpha_num + d_num, denum))
# # # # ]
43 : 0 : return 1;
44 : :
45 : : return 0;
46 : : }
47 : :
48 : : static inline void
49 : : find_exact_solution_left(uint32_t p_a, uint32_t q_a, uint32_t p_b, uint32_t q_b,
50 : : uint32_t alpha_num, uint32_t d_num, uint32_t denum, uint32_t *p, uint32_t *q)
51 : : {
52 : 0 : uint32_t k_num = denum * p_b - (alpha_num + d_num) * q_b;
53 : 0 : uint32_t k_denum = (alpha_num + d_num) * q_a - denum * p_a;
54 : 0 : uint32_t k = (k_num / k_denum) + 1;
55 : :
56 : 0 : *p = p_b + k * p_a;
57 : 0 : *q = q_b + k * q_a;
58 : : }
59 : :
60 : : static inline void
61 : : find_exact_solution_right(uint32_t p_a, uint32_t q_a, uint32_t p_b, uint32_t q_b,
62 : : uint32_t alpha_num, uint32_t d_num, uint32_t denum, uint32_t *p, uint32_t *q)
63 : : {
64 : 0 : uint32_t k_num = - denum * p_b + (alpha_num - d_num) * q_b;
65 : 0 : uint32_t k_denum = - (alpha_num - d_num) * q_a + denum * p_a;
66 : 0 : uint32_t k = (k_num / k_denum) + 1;
67 : :
68 : 0 : *p = p_b + k * p_a;
69 : 0 : *q = q_b + k * q_a;
70 : : }
71 : :
72 : : static int
73 : 0 : find_best_rational_approximation(uint32_t alpha_num, uint32_t d_num, uint32_t denum, uint32_t *p, uint32_t *q)
74 : : {
75 : : uint32_t p_a, q_a, p_b, q_b;
76 : :
77 : : /* check assumptions on the inputs */
78 [ # # # # : 0 : if (!((0 < d_num) && (d_num < alpha_num) && (alpha_num < denum) && (d_num + alpha_num < denum))) {
# # ]
79 : : return -1;
80 : : }
81 : :
82 : : /* set initial bounds for the search */
83 : : p_a = 0;
84 : : q_a = 1;
85 : : p_b = 1;
86 : : q_b = 1;
87 : :
88 : : while (1) {
89 : : uint32_t new_p_a, new_q_a, new_p_b, new_q_b;
90 : : uint32_t x_num, x_denum, x;
91 : : int aa, bb;
92 : :
93 : : /* compute the number of steps to the left */
94 : 0 : x_num = denum * p_b - alpha_num * q_b;
95 : 0 : x_denum = - denum * p_a + alpha_num * q_a;
96 : 0 : x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */
97 : :
98 : : /* check whether we have a valid approximation */
99 [ # # ]: 0 : aa = matches(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum);
100 [ # # ]: 0 : bb = matches(p_b + (x-1) * p_a, q_b + (x - 1) * q_a, alpha_num, d_num, denum);
101 [ # # ]: 0 : if (aa || bb) {
102 : : find_exact_solution_left(p_a, q_a, p_b, q_b, alpha_num, d_num, denum, p, q);
103 : 0 : return 0;
104 : : }
105 : :
106 : : /* update the interval */
107 : : new_p_a = p_b + (x - 1) * p_a ;
108 : : new_q_a = q_b + (x - 1) * q_a;
109 : : new_p_b = p_b + x * p_a ;
110 : : new_q_b = q_b + x * q_a;
111 : :
112 : : p_a = new_p_a ;
113 : : q_a = new_q_a;
114 : : p_b = new_p_b ;
115 : : q_b = new_q_b;
116 : :
117 : : /* compute the number of steps to the right */
118 : 0 : x_num = alpha_num * q_b - denum * p_b;
119 : 0 : x_denum = - alpha_num * q_a + denum * p_a;
120 : 0 : x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */
121 : :
122 : : /* check whether we have a valid approximation */
123 [ # # ]: 0 : aa = matches(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum);
124 [ # # ]: 0 : bb = matches(p_b + (x - 1) * p_a, q_b + (x - 1) * q_a, alpha_num, d_num, denum);
125 [ # # ]: 0 : if (aa || bb) {
126 : : find_exact_solution_right(p_a, q_a, p_b, q_b, alpha_num, d_num, denum, p, q);
127 : 0 : return 0;
128 : : }
129 : :
130 : : /* update the interval */
131 : : new_p_a = p_b + (x - 1) * p_a;
132 : : new_q_a = q_b + (x - 1) * q_a;
133 : : new_p_b = p_b + x * p_a;
134 : : new_q_b = q_b + x * q_a;
135 : :
136 : : p_a = new_p_a;
137 : : q_a = new_q_a;
138 : : p_b = new_p_b;
139 : : q_b = new_q_b;
140 : : }
141 : : }
142 : :
143 : : RTE_EXPORT_SYMBOL(rte_approx)
144 : 0 : int rte_approx(double alpha, double d, uint32_t *p, uint32_t *q)
145 : : {
146 : : uint32_t alpha_num, d_num, denum;
147 : :
148 : : /* Check input arguments */
149 [ # # # # : 0 : if (!((0.0 < d) && (d < alpha) && (alpha < 1.0))) {
# # ]
150 : : return -1;
151 : : }
152 : :
153 [ # # ]: 0 : if ((p == NULL) || (q == NULL)) {
154 : : return -2;
155 : : }
156 : :
157 : : /* Compute alpha_num, d_num and denum */
158 : : denum = 1;
159 [ # # ]: 0 : while (d < 1) {
160 : 0 : alpha *= 10;
161 : 0 : d *= 10;
162 : 0 : denum *= 10;
163 : : }
164 : 0 : alpha_num = (uint32_t) alpha;
165 : 0 : d_num = (uint32_t) d;
166 : :
167 : : /* Perform approximation */
168 : 0 : return find_best_rational_approximation(alpha_num, d_num, denum, p, q);
169 : : }
170 : :
171 : : /* fraction comparison (64 bit version): compare (a/b) and (c/d) */
172 : : static inline uint64_t
173 : : less_64(uint64_t a, uint64_t b, uint64_t c, uint64_t d)
174 : : {
175 : 2 : return a*d < b*c;
176 : : }
177 : :
178 : : static inline uint64_t
179 : : less_or_equal_64(uint64_t a, uint64_t b, uint64_t c, uint64_t d)
180 : : {
181 : 2 : return a*d <= b*c;
182 : : }
183 : :
184 : : /* check whether a/b is a valid approximation (64 bit version) */
185 : : static inline uint64_t
186 : : matches_64(uint64_t a, uint64_t b,
187 : : uint64_t alpha_num, uint64_t d_num, uint64_t denum)
188 : : {
189 : 2 : if (less_or_equal_64(a, b, alpha_num - d_num, denum))
190 : : return 0;
191 : :
192 [ + - + - : 2 : if (less_64(a, b, alpha_num + d_num, denum))
- - - - ]
193 : 2 : return 1;
194 : :
195 : : return 0;
196 : : }
197 : :
198 : : static inline void
199 : : find_exact_solution_left_64(uint64_t p_a, uint64_t q_a, uint64_t p_b, uint64_t q_b,
200 : : uint64_t alpha_num, uint64_t d_num, uint64_t denum, uint64_t *p, uint64_t *q)
201 : : {
202 : 1 : uint64_t k_num = denum * p_b - (alpha_num + d_num) * q_b;
203 : 1 : uint64_t k_denum = (alpha_num + d_num) * q_a - denum * p_a;
204 : 1 : uint64_t k = (k_num / k_denum) + 1;
205 : :
206 : 1 : *p = p_b + k * p_a;
207 : 1 : *q = q_b + k * q_a;
208 : : }
209 : :
210 : : static inline void
211 : : find_exact_solution_right_64(uint64_t p_a, uint64_t q_a, uint64_t p_b, uint64_t q_b,
212 : : uint64_t alpha_num, uint64_t d_num, uint64_t denum, uint64_t *p, uint64_t *q)
213 : : {
214 : 0 : uint64_t k_num = -denum * p_b + (alpha_num - d_num) * q_b;
215 : 0 : uint64_t k_denum = -(alpha_num - d_num) * q_a + denum * p_a;
216 : 0 : uint64_t k = (k_num / k_denum) + 1;
217 : :
218 : 0 : *p = p_b + k * p_a;
219 : 0 : *q = q_b + k * q_a;
220 : : }
221 : :
222 : : static int
223 : 1 : find_best_rational_approximation_64(uint64_t alpha_num, uint64_t d_num,
224 : : uint64_t denum, uint64_t *p, uint64_t *q)
225 : : {
226 : : uint64_t p_a, q_a, p_b, q_b;
227 : :
228 : : /* check assumptions on the inputs */
229 [ + - + - ]: 1 : if (!((d_num > 0) && (d_num < alpha_num) &&
230 [ + - ]: 1 : (alpha_num < denum) && (d_num + alpha_num < denum))) {
231 : : return -1;
232 : : }
233 : :
234 : : /* set initial bounds for the search */
235 : : p_a = 0;
236 : : q_a = 1;
237 : : p_b = 1;
238 : : q_b = 1;
239 : :
240 : : while (1) {
241 : : uint64_t new_p_a, new_q_a, new_p_b, new_q_b;
242 : : uint64_t x_num, x_denum, x;
243 : : int aa, bb;
244 : :
245 : : /* compute the number of steps to the left */
246 : 1 : x_num = denum * p_b - alpha_num * q_b;
247 : 1 : x_denum = -denum * p_a + alpha_num * q_a;
248 : 1 : x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */
249 : :
250 : : /* check whether we have a valid approximation */
251 [ + - ]: 1 : aa = matches_64(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum);
252 [ + - ]: 1 : bb = matches_64(p_b + (x-1) * p_a, q_b + (x - 1) * q_a,
253 : : alpha_num, d_num, denum);
254 [ + - ]: 1 : if (aa || bb) {
255 : : find_exact_solution_left_64(p_a, q_a, p_b, q_b,
256 : : alpha_num, d_num, denum, p, q);
257 : 1 : return 0;
258 : : }
259 : :
260 : : /* update the interval */
261 : : new_p_a = p_b + (x - 1) * p_a;
262 : : new_q_a = q_b + (x - 1) * q_a;
263 : : new_p_b = p_b + x * p_a;
264 : : new_q_b = q_b + x * q_a;
265 : :
266 : : p_a = new_p_a;
267 : : q_a = new_q_a;
268 : : p_b = new_p_b;
269 : : q_b = new_q_b;
270 : :
271 : : /* compute the number of steps to the right */
272 : 0 : x_num = alpha_num * q_b - denum * p_b;
273 : 0 : x_denum = -alpha_num * q_a + denum * p_a;
274 : 0 : x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */
275 : :
276 : : /* check whether we have a valid approximation */
277 [ # # ]: 0 : aa = matches_64(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum);
278 [ # # ]: 0 : bb = matches_64(p_b + (x - 1) * p_a, q_b + (x - 1) * q_a,
279 : : alpha_num, d_num, denum);
280 [ # # ]: 0 : if (aa || bb) {
281 : : find_exact_solution_right_64(p_a, q_a, p_b, q_b,
282 : : alpha_num, d_num, denum, p, q);
283 : 0 : return 0;
284 : : }
285 : :
286 : : /* update the interval */
287 : : new_p_a = p_b + (x - 1) * p_a;
288 : : new_q_a = q_b + (x - 1) * q_a;
289 : : new_p_b = p_b + x * p_a;
290 : : new_q_b = q_b + x * q_a;
291 : :
292 : : p_a = new_p_a;
293 : : q_a = new_q_a;
294 : : p_b = new_p_b;
295 : : q_b = new_q_b;
296 : : }
297 : : }
298 : :
299 : 1 : int rte_approx_64(double alpha, double d, uint64_t *p, uint64_t *q)
300 : : {
301 : : uint64_t alpha_num, d_num, denum;
302 : :
303 : : /* Check input arguments */
304 [ + - + - : 1 : if (!((0.0 < d) && (d < alpha) && (alpha < 1.0)))
+ - ]
305 : : return -1;
306 : :
307 [ + - ]: 1 : if ((p == NULL) || (q == NULL))
308 : : return -2;
309 : :
310 : : /* Compute alpha_num, d_num and denum */
311 : : denum = 1;
312 [ + + ]: 8 : while (d < 1) {
313 : 7 : alpha *= 10;
314 : 7 : d *= 10;
315 : 7 : denum *= 10;
316 : : }
317 : 1 : alpha_num = (uint64_t) alpha;
318 : 1 : d_num = (uint64_t) d;
319 : :
320 : : /* Perform approximation */
321 : 1 : return find_best_rational_approximation_64(alpha_num, d_num, denum, p, q);
322 : : }
|
