| Vinod Polimera | 7528311 | 2022-07-20 17:25:44 +0530 | [diff] [blame] | 1 | /* SPDX-License-Identifier: GPL-2.0-only */ |
| 2 | /* |
| 3 | * Helper functions for rational numbers. |
| 4 | * |
| 5 | * Copyright (C) 2009 emlix GmbH, Oskar Schirmer <oskar@scara.com> |
| 6 | * Copyright (C) 2019 Trent Piepho <tpiepho@gmail.com> |
| 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 | */ |
| 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 | } |