lib/math/rational.c

0001 // SPDX-License-Identifier: GPL-2.0
0002 /*
0003  * rational fractions
0004  *
0005  * Copyright (C) 2009 emlix GmbH, Oskar Schirmer <oskar@scara.com>
0006  * Copyright (C) 2019 Trent Piepho <tpiepho@gmail.com>
0007  *
0008  * helper functions when coping with rational numbers
0009  */
0010
0011 #include <linux/rational.h>
0012 #include <linux/compiler.h>
0013 #include <linux/export.h>
0014 #include <linux/minmax.h>
0015 #include <linux/limits.h>
0016 #include <linux/module.h>
0017
0018 /*
0019  * calculate best rational approximation for a given fraction
0020  * taking into account restricted register size, e.g. to find
0021  * appropriate values for a pll with 5 bit denominator and
0022  * 8 bit numerator register fields, trying to set up with a
0023  * frequency ratio of 3.1415, one would say:
0024  *
0025  * rational_best_approximation(31415, 10000,
0026  *      (1 << 8) - 1, (1 << 5) - 1, &n, &d);
0027  *
0028  * you may look at given_numerator as a fixed point number,
0029  * with the fractional part size described in given_denominator.
0030  *
0031  * for theoretical background, see:
0032  * https://en.wikipedia.org/wiki/Continued_fraction
0033  */
0034
0035 void rational_best_approximation(
0036     unsigned long given_numerator, unsigned long given_denominator,
0037     unsigned long max_numerator, unsigned long max_denominator,
0038     unsigned long *best_numerator, unsigned long *best_denominator)
0039 {
0040     /* n/d is the starting rational, which is continually
0041      * decreased each iteration using the Euclidean algorithm.
0042      *
0043      * dp is the value of d from the prior iteration.
0044      *
0045      * n2/d2, n1/d1, and n0/d0 are our successively more accurate
0046      * approximations of the rational.  They are, respectively,
0047      * the current, previous, and two prior iterations of it.
0048      *
0049      * a is current term of the continued fraction.
0050      */
0051     unsigned long n, d, n0, d0, n1, d1, n2, d2;
0052     n = given_numerator;
0053     d = given_denominator;
0054     n0 = d1 = 0;
0055     n1 = d0 = 1;
0056
0057     for (;;) {
0058         unsigned long dp, a;
0059
0060         if (d == 0)
0061             break;
0062         /* Find next term in continued fraction, 'a', via
0063          * Euclidean algorithm.
0064          */
0065         dp = d;
0066         a = n / d;
0067         d = n % d;
0068         n = dp;
0069
0070         /* Calculate the current rational approximation (aka
0071          * convergent), n2/d2, using the term just found and
0072          * the two prior approximations.
0073          */
0074         n2 = n0 + a * n1;
0075         d2 = d0 + a * d1;
0076
0077         /* If the current convergent exceeds the maxes, then
0078          * return either the previous convergent or the
0079          * largest semi-convergent, the final term of which is
0080          * found below as 't'.
0081          */
0082         if ((n2 > max_numerator) || (d2 > max_denominator)) {
0083             unsigned long t = ULONG_MAX;
0084
0085             if (d1)
0086                 t = (max_denominator - d0) / d1;
0087             if (n1)
0088                 t = min(t, (max_numerator - n0) / n1);
0089
0090             /* This tests if the semi-convergent is closer than the previous
0091              * convergent.  If d1 is zero there is no previous convergent as this
0092              * is the 1st iteration, so always choose the semi-convergent.
0093              */
0094             if (!d1 || 2u * t > a || (2u * t == a && d0 * dp > d1 * d)) {
0095                 n1 = n0 + t * n1;
0096                 d1 = d0 + t * d1;
0097             }
0098             break;
0099         }
0100         n0 = n1;
0101         n1 = n2;
0102         d0 = d1;
0103         d1 = d2;
0104     }
0105     *best_numerator = n1;
0106     *best_denominator = d1;
0107 }
0108
0109 EXPORT_SYMBOL(rational_best_approximation);
0110
0111 MODULE_LICENSE("GPL v2");