Branch data     Line data    Source code

       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/kernel.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                 :            :  * http://en.wikipedia.org/wiki/Continued_fraction
      31                 :            :  */
      32                 :            : 
      33                 :          0 : 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                 :          0 :         unsigned long n, d, n0, d0, n1, d1, n2, d2;
      50                 :          0 :         n = given_numerator;
      51                 :          0 :         d = given_denominator;
      52                 :          0 :         n0 = d1 = 0;
      53                 :          0 :         n1 = d0 = 1;
      54                 :            : 
      55                 :          0 :         for (;;) {
      56                 :          0 :                 unsigned long dp, a;
      57                 :            : 
      58         [ #  # ]:          0 :                 if (d == 0)
      59                 :            :                         break;
      60                 :            :                 /* Find next term in continued fraction, 'a', via
      61                 :            :                  * Euclidean algorithm.
      62                 :            :                  */
      63                 :          0 :                 dp = d;
      64                 :          0 :                 a = n / d;
      65                 :          0 :                 d = n % d;
      66                 :          0 :                 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                 :          0 :                 n2 = n0 + a * n1;
      73                 :          0 :                 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         [ #  # ]:          0 :                 if ((n2 > max_numerator) || (d2 > max_denominator)) {
      81                 :          0 :                         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   [ #  #  #  #  :          0 :                         if (2u * t > a || (2u * t == a && d0 * dp > d1 * d)) {
                   #  # ]
      88                 :          0 :                                 n1 = n0 + t * n1;
      89                 :          0 :                                 d1 = d0 + t * d1;
      90                 :            :                         }
      91                 :            :                         break;
      92                 :            :                 }
      93                 :            :                 n0 = n1;
      94                 :            :                 n1 = n2;
      95                 :            :                 d0 = d1;
      96                 :            :                 d1 = d2;
      97                 :            :         }
      98                 :          0 :         *best_numerator = n1;
      99                 :          0 :         *best_denominator = d1;
     100                 :          0 : }
     101                 :            : 
     102                 :            : EXPORT_SYMBOL(rational_best_approximation);