1 // SPDX-License-Identifier: GPL-2.0 2 /* 3 * rational fractions 4 * 5 * Copyright (C) 2009 emlix GmbH, Oskar Schirmer <oskar@scara.com> 6 * Copyright (C) 2019 Trent Piepho <tpiepho@gmail.com> 7 * 8 * helper functions when coping with rational numbers 9 */ 10 11 #include <linux/rational.h> 12 #include <linux/compiler.h> 13 #include <linux/export.h> 14 #include <linux/minmax.h> 15 16 /* 17 * calculate best rational approximation for a given fraction 18 * taking into account restricted register size, e.g. to find 19 * appropriate values for a pll with 5 bit denominator and 20 * 8 bit numerator register fields, trying to set up with a 21 * frequency ratio of 3.1415, one would say: 22 * 23 * rational_best_approximation(31415, 10000, 24 * (1 << 8) - 1, (1 << 5) - 1, &n, &d); 25 * 26 * you may look at given_numerator as a fixed point number, 27 * with the fractional part size described in given_denominator. 28 * 29 * for theoretical background, see: 30 * https://en.wikipedia.org/wiki/Continued_fraction 31 */ 32 33 void rational_best_approximation( 34 unsigned long given_numerator, unsigned long given_denominator, 35 unsigned long max_numerator, unsigned long max_denominator, 36 unsigned long *best_numerator, unsigned long *best_denominator) 37 { 38 /* n/d is the starting rational, which is continually 39 * decreased each iteration using the Euclidean algorithm. 40 * 41 * dp is the value of d from the prior iteration. 42 * 43 * n2/d2, n1/d1, and n0/d0 are our successively more accurate 44 * approximations of the rational. They are, respectively, 45 * the current, previous, and two prior iterations of it. 46 * 47 * a is current term of the continued fraction. 48 */ 49 unsigned long n, d, n0, d0, n1, d1, n2, d2; 50 n = given_numerator; 51 d = given_denominator; 52 n0 = d1 = 0; 53 n1 = d0 = 1; 54 55 for (;;) { 56 unsigned long dp, a; 57 58 if (d == 0) 59 break; 60 /* Find next term in continued fraction, 'a', via 61 * Euclidean algorithm. 62 */ 63 dp = d; 64 a = n / d; 65 d = n % d; 66 n = dp; 67 68 /* Calculate the current rational approximation (aka 69 * convergent), n2/d2, using the term just found and 70 * the two prior approximations. 71 */ 72 n2 = n0 + a * n1; 73 d2 = d0 + a * d1; 74 75 /* If the current convergent exceeds the maxes, then 76 * return either the previous convergent or the 77 * largest semi-convergent, the final term of which is 78 * found below as 't'. 79 */ 80 if ((n2 > max_numerator) || (d2 > max_denominator)) { 81 unsigned long t = min((max_numerator - n0) / n1, 82 (max_denominator - d0) / d1); 83 84 /* This tests if the semi-convergent is closer 85 * than the previous convergent. 86 */ 87 if (2u * t > a || (2u * t == a && d0 * dp > d1 * d)) { 88 n1 = n0 + t * n1; 89 d1 = d0 + t * d1; 90 } 91 break; 92 } 93 n0 = n1; 94 n1 = n2; 95 d0 = d1; 96 d1 = d2; 97 } 98 *best_numerator = n1; 99 *best_denominator = d1; 100 } 101 102 EXPORT_SYMBOL(rational_best_approximation); 103