1 /*- 2 * BSD LICENSE 3 * 4 * Copyright(c) 2010-2014 Intel Corporation. All rights reserved. 5 * All rights reserved. 6 * 7 * Redistribution and use in source and binary forms, with or without 8 * modification, are permitted provided that the following conditions 9 * are met: 10 * 11 * * Redistributions of source code must retain the above copyright 12 * notice, this list of conditions and the following disclaimer. 13 * * Redistributions in binary form must reproduce the above copyright 14 * notice, this list of conditions and the following disclaimer in 15 * the documentation and/or other materials provided with the 16 * distribution. 17 * * Neither the name of Intel Corporation nor the names of its 18 * contributors may be used to endorse or promote products derived 19 * from this software without specific prior written permission. 20 * 21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 24 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 25 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 26 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 27 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 28 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 29 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 30 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 31 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 32 */ 33 34 #include <stdlib.h> 35 36 #include "rte_approx.h" 37 38 /* 39 * Based on paper "Approximating Rational Numbers by Fractions" by Michal 40 * Forisek [email protected] 41 * 42 * Given a rational number alpha with 0 < alpha < 1 and a precision d, the goal 43 * is to find positive integers p, q such that alpha - d < p/q < alpha + d, and 44 * q is minimal. 45 * 46 * http://people.ksp.sk/~misof/publications/2007approx.pdf 47 */ 48 49 /* fraction comparison: compare (a/b) and (c/d) */ 50 static inline uint32_t 51 less(uint32_t a, uint32_t b, uint32_t c, uint32_t d) 52 { 53 return a*d < b*c; 54 } 55 56 static inline uint32_t 57 less_or_equal(uint32_t a, uint32_t b, uint32_t c, uint32_t d) 58 { 59 return a*d <= b*c; 60 } 61 62 /* check whether a/b is a valid approximation */ 63 static inline uint32_t 64 matches(uint32_t a, uint32_t b, 65 uint32_t alpha_num, uint32_t d_num, uint32_t denum) 66 { 67 if (less_or_equal(a, b, alpha_num - d_num, denum)) 68 return 0; 69 70 if (less(a ,b, alpha_num + d_num, denum)) 71 return 1; 72 73 return 0; 74 } 75 76 static inline void 77 find_exact_solution_left(uint32_t p_a, uint32_t q_a, uint32_t p_b, uint32_t q_b, 78 uint32_t alpha_num, uint32_t d_num, uint32_t denum, uint32_t *p, uint32_t *q) 79 { 80 uint32_t k_num = denum * p_b - (alpha_num + d_num) * q_b; 81 uint32_t k_denum = (alpha_num + d_num) * q_a - denum * p_a; 82 uint32_t k = (k_num / k_denum) + 1; 83 84 *p = p_b + k * p_a; 85 *q = q_b + k * q_a; 86 } 87 88 static inline void 89 find_exact_solution_right(uint32_t p_a, uint32_t q_a, uint32_t p_b, uint32_t q_b, 90 uint32_t alpha_num, uint32_t d_num, uint32_t denum, uint32_t *p, uint32_t *q) 91 { 92 uint32_t k_num = - denum * p_b + (alpha_num - d_num) * q_b; 93 uint32_t k_denum = - (alpha_num - d_num) * q_a + denum * p_a; 94 uint32_t k = (k_num / k_denum) + 1; 95 96 *p = p_b + k * p_a; 97 *q = q_b + k * q_a; 98 } 99 100 static int 101 find_best_rational_approximation(uint32_t alpha_num, uint32_t d_num, uint32_t denum, uint32_t *p, uint32_t *q) 102 { 103 uint32_t p_a, q_a, p_b, q_b; 104 105 /* check assumptions on the inputs */ 106 if (!((0 < d_num) && (d_num < alpha_num) && (alpha_num < denum) && (d_num + alpha_num < denum))) { 107 return -1; 108 } 109 110 /* set initial bounds for the search */ 111 p_a = 0; 112 q_a = 1; 113 p_b = 1; 114 q_b = 1; 115 116 while (1) { 117 uint32_t new_p_a, new_q_a, new_p_b, new_q_b; 118 uint32_t x_num, x_denum, x; 119 int aa, bb; 120 121 /* compute the number of steps to the left */ 122 x_num = denum * p_b - alpha_num * q_b; 123 x_denum = - denum * p_a + alpha_num * q_a; 124 x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */ 125 126 /* check whether we have a valid approximation */ 127 aa = matches(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum); 128 bb = matches(p_b + (x-1) * p_a, q_b + (x - 1) * q_a, alpha_num, d_num, denum); 129 if (aa || bb) { 130 find_exact_solution_left(p_a, q_a, p_b, q_b, alpha_num, d_num, denum, p, q); 131 return 0; 132 } 133 134 /* update the interval */ 135 new_p_a = p_b + (x - 1) * p_a ; 136 new_q_a = q_b + (x - 1) * q_a; 137 new_p_b = p_b + x * p_a ; 138 new_q_b = q_b + x * q_a; 139 140 p_a = new_p_a ; 141 q_a = new_q_a; 142 p_b = new_p_b ; 143 q_b = new_q_b; 144 145 /* compute the number of steps to the right */ 146 x_num = alpha_num * q_b - denum * p_b; 147 x_denum = - alpha_num * q_a + denum * p_a; 148 x = (x_num + x_denum - 1) / x_denum; /* x = ceil(x_num / x_denum) */ 149 150 /* check whether we have a valid approximation */ 151 aa = matches(p_b + x * p_a, q_b + x * q_a, alpha_num, d_num, denum); 152 bb = matches(p_b + (x - 1) * p_a, q_b + (x - 1) * q_a, alpha_num, d_num, denum); 153 if (aa || bb) { 154 find_exact_solution_right(p_a, q_a, p_b, q_b, alpha_num, d_num, denum, p, q); 155 return 0; 156 } 157 158 /* update the interval */ 159 new_p_a = p_b + (x - 1) * p_a; 160 new_q_a = q_b + (x - 1) * q_a; 161 new_p_b = p_b + x * p_a; 162 new_q_b = q_b + x * q_a; 163 164 p_a = new_p_a; 165 q_a = new_q_a; 166 p_b = new_p_b; 167 q_b = new_q_b; 168 } 169 } 170 171 int rte_approx(double alpha, double d, uint32_t *p, uint32_t *q) 172 { 173 uint32_t alpha_num, d_num, denum; 174 175 /* Check input arguments */ 176 if (!((0.0 < d) && (d < alpha) && (alpha < 1.0))) { 177 return -1; 178 } 179 180 if ((p == NULL) || (q == NULL)) { 181 return -2; 182 } 183 184 /* Compute alpha_num, d_num and denum */ 185 denum = 1; 186 while (d < 1) { 187 alpha *= 10; 188 d *= 10; 189 denum *= 10; 190 } 191 alpha_num = (uint32_t) alpha; 192 d_num = (uint32_t) d; 193 194 /* Perform approximation */ 195 return find_best_rational_approximation(alpha_num, d_num, denum, p, q); 196 } 197