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);
|