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 [email protected] 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