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 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 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 if (less_or_equal(a, b, alpha_num - d_num, denum)) 39 return 0; 40 41 if (less(a ,b, alpha_num + d_num, denum)) 42 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 uint32_t k_num = denum * p_b - (alpha_num + d_num) * q_b; 52 uint32_t k_denum = (alpha_num + d_num) * q_a - denum * p_a; 53 uint32_t k = (k_num / k_denum) + 1; 54 55 *p = p_b + k * p_a; 56 *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 uint32_t k_num = - denum * p_b + (alpha_num - d_num) * q_b; 64 uint32_t k_denum = - (alpha_num - d_num) * q_a + denum * p_a; 65 uint32_t k = (k_num / k_denum) + 1; 66 67 *p = p_b + k * p_a; 68 *q = q_b + k * q_a; 69 } 70 71 static int 72 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 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 x_num = denum * p_b - alpha_num * q_b; 94 x_denum = - denum * p_a + alpha_num * q_a; 95 x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */ 96 97 /* check whether we have a valid approximation */ 98 aa = matches(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum); 99 bb = matches(p_b + (x-1) * p_a, q_b + (x - 1) * q_a, alpha_num, d_num, denum); 100 if (aa || bb) { 101 find_exact_solution_left(p_a, q_a, p_b, q_b, alpha_num, d_num, denum, p, q); 102 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 x_num = alpha_num * q_b - denum * p_b; 118 x_denum = - alpha_num * q_a + denum * p_a; 119 x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */ 120 121 /* check whether we have a valid approximation */ 122 aa = matches(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum); 123 bb = matches(p_b + (x - 1) * p_a, q_b + (x - 1) * q_a, alpha_num, d_num, denum); 124 if (aa || bb) { 125 find_exact_solution_right(p_a, q_a, p_b, q_b, alpha_num, d_num, denum, p, q); 126 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 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 if (!((0.0 < d) && (d < alpha) && (alpha < 1.0))) { 148 return -1; 149 } 150 151 if ((p == NULL) || (q == NULL)) { 152 return -2; 153 } 154 155 /* Compute alpha_num, d_num and denum */ 156 denum = 1; 157 while (d < 1) { 158 alpha *= 10; 159 d *= 10; 160 denum *= 10; 161 } 162 alpha_num = (uint32_t) alpha; 163 d_num = (uint32_t) d; 164 165 /* Perform approximation */ 166 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 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 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 if (less_or_equal_64(a, b, alpha_num - d_num, denum)) 188 return 0; 189 190 if (less_64(a, b, alpha_num + d_num, denum)) 191 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 uint64_t k_num = denum * p_b - (alpha_num + d_num) * q_b; 201 uint64_t k_denum = (alpha_num + d_num) * q_a - denum * p_a; 202 uint64_t k = (k_num / k_denum) + 1; 203 204 *p = p_b + k * p_a; 205 *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 uint64_t k_num = -denum * p_b + (alpha_num - d_num) * q_b; 213 uint64_t k_denum = -(alpha_num - d_num) * q_a + denum * p_a; 214 uint64_t k = (k_num / k_denum) + 1; 215 216 *p = p_b + k * p_a; 217 *q = q_b + k * q_a; 218 } 219 220 static int 221 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 if (!((d_num > 0) && (d_num < alpha_num) && 228 (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 x_num = denum * p_b - alpha_num * q_b; 245 x_denum = -denum * p_a + alpha_num * q_a; 246 x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */ 247 248 /* check whether we have a valid approximation */ 249 aa = matches_64(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum); 250 bb = matches_64(p_b + (x-1) * p_a, q_b + (x - 1) * q_a, 251 alpha_num, d_num, denum); 252 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 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 x_num = alpha_num * q_b - denum * p_b; 271 x_denum = -alpha_num * q_a + denum * p_a; 272 x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */ 273 274 /* check whether we have a valid approximation */ 275 aa = matches_64(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum); 276 bb = matches_64(p_b + (x - 1) * p_a, q_b + (x - 1) * q_a, 277 alpha_num, d_num, denum); 278 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 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 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 if (!((0.0 < d) && (d < alpha) && (alpha < 1.0))) 303 return -1; 304 305 if ((p == NULL) || (q == NULL)) 306 return -2; 307 308 /* Compute alpha_num, d_num and denum */ 309 denum = 1; 310 while (d < 1) { 311 alpha *= 10; 312 d *= 10; 313 denum *= 10; 314 } 315 alpha_num = (uint64_t) alpha; 316 d_num = (uint64_t) d; 317 318 /* Perform approximation */ 319 return find_best_rational_approximation_64(alpha_num, d_num, denum, p, q); 320 } 321