1 /* SPDX-License-Identifier: GPL-2.0-only */
2 /*
3 * Helper functions for rational numbers.
4 *
5 * Copyright (C) 2009 emlix GmbH, Oskar Schirmer <[email protected]>
6 * Copyright (C) 2019 Trent Piepho <[email protected]>
7 */
8
9 #include <commonlib/helpers.h>
10 #include <commonlib/rational.h>
11 #include <limits.h>
12
13 /*
14 * For theoretical background, see:
15 * https://en.wikipedia.org/wiki/Continued_fraction
16 */
rational_best_approximation(unsigned long numerator,unsigned long denominator,unsigned long max_numerator,unsigned long max_denominator,unsigned long * best_numerator,unsigned long * best_denominator)17 void rational_best_approximation(
18 unsigned long numerator, unsigned long denominator,
19 unsigned long max_numerator, unsigned long max_denominator,
20 unsigned long *best_numerator, unsigned long *best_denominator)
21 {
22 /*
23 * n/d is the starting rational, where both n and d will
24 * decrease in each iteration using the Euclidean algorithm.
25 *
26 * dp is the value of d from the prior iteration.
27 *
28 * n2/d2, n1/d1, and n0/d0 are our successively more accurate
29 * approximations of the rational. They are, respectively,
30 * the current, previous, and two prior iterations of it.
31 *
32 * a is current term of the continued fraction.
33 */
34 unsigned long n, d, n0, d0, n1, d1, n2, d2;
35 n = numerator;
36 d = denominator;
37 n0 = d1 = 0;
38 n1 = d0 = 1;
39
40 for (;;) {
41 unsigned long dp, a;
42
43 if (d == 0)
44 break;
45 /*
46 * Find next term in continued fraction, 'a', via
47 * Euclidean algorithm.
48 */
49 dp = d;
50 a = n / d;
51 d = n % d;
52 n = dp;
53
54 /*
55 * Calculate the current rational approximation (aka
56 * convergent), n2/d2, using the term just found and
57 * the two prior approximations.
58 */
59 n2 = n0 + a * n1;
60 d2 = d0 + a * d1;
61
62 /*
63 * If the current convergent exceeds the maximum, then
64 * return either the previous convergent or the
65 * largest semi-convergent, the final term of which is
66 * found below as 't'.
67 */
68 if ((n2 > max_numerator) || (d2 > max_denominator)) {
69 unsigned long t = ULONG_MAX;
70
71 if (d1)
72 t = (max_denominator - d0) / d1;
73 if (n1)
74 t = MIN(t, (max_numerator - n0) / n1);
75
76 /*
77 * This tests if the semi-convergent is closer than the previous
78 * convergent. If d1 is zero there is no previous convergent as
79 * this is the 1st iteration, so always choose the semi-convergent.
80 */
81 if (!d1 || 2u * t > a || (2u * t == a && d0 * dp > d1 * d)) {
82 n1 = n0 + t * n1;
83 d1 = d0 + t * d1;
84 }
85 break;
86 }
87 n0 = n1;
88 n1 = n2;
89 d0 = d1;
90 d1 = d2;
91 }
92
93 *best_numerator = n1;
94 *best_denominator = d1;
95 }
96