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