linux-stable/lib/math/rational.c
Trent Piepho 65a0d3c146 lib/math/rational.c: fix divide by zero
If the input is out of the range of the allowed values, either larger than
the largest value or closer to zero than the smallest non-zero allowed
value, then a division by zero would occur.

In the case of input too large, the division by zero will occur on the
first iteration.  The best result (largest allowed value) will be found by
always choosing the semi-convergent and excluding the denominator based
limit when finding it.

In the case of the input too small, the division by zero will occur on the
second iteration.  The numerator based semi-convergent should not be
calculated to avoid the division by zero.  But the semi-convergent vs
previous convergent test is still needed, which effectively chooses
between 0 (the previous convergent) vs the smallest allowed fraction (best
semi-convergent) as the result.

Link: https://lkml.kernel.org/r/20210525144250.214670-1-tpiepho@gmail.com
Fixes: 323dd2c3ed ("lib/math/rational.c: fix possible incorrect result from rational fractions helper")
Signed-off-by: Trent Piepho <tpiepho@gmail.com>
Reported-by: Yiyuan Guo <yguoaz@gmail.com>
Reviewed-by: Andy Shevchenko <andriy.shevchenko@linux.intel.com>
Cc: Oskar Schirmer <oskar@scara.com>
Cc: Daniel Latypov <dlatypov@google.com>
Signed-off-by: Andrew Morton <akpm@linux-foundation.org>
Signed-off-by: Linus Torvalds <torvalds@linux-foundation.org>
2021-07-01 11:06:05 -07:00

108 lines
2.9 KiB
C

// SPDX-License-Identifier: GPL-2.0
/*
* rational fractions
*
* Copyright (C) 2009 emlix GmbH, Oskar Schirmer <oskar@scara.com>
* Copyright (C) 2019 Trent Piepho <tpiepho@gmail.com>
*
* helper functions when coping with rational numbers
*/
#include <linux/rational.h>
#include <linux/compiler.h>
#include <linux/export.h>
#include <linux/minmax.h>
#include <linux/limits.h>
/*
* calculate best rational approximation for a given fraction
* taking into account restricted register size, e.g. to find
* appropriate values for a pll with 5 bit denominator and
* 8 bit numerator register fields, trying to set up with a
* frequency ratio of 3.1415, one would say:
*
* rational_best_approximation(31415, 10000,
* (1 << 8) - 1, (1 << 5) - 1, &n, &d);
*
* you may look at given_numerator as a fixed point number,
* with the fractional part size described in given_denominator.
*
* for theoretical background, see:
* https://en.wikipedia.org/wiki/Continued_fraction
*/
void rational_best_approximation(
unsigned long given_numerator, unsigned long given_denominator,
unsigned long max_numerator, unsigned long max_denominator,
unsigned long *best_numerator, unsigned long *best_denominator)
{
/* n/d is the starting rational, which is continually
* decreased each iteration using the Euclidean algorithm.
*
* dp is the value of d from the prior iteration.
*
* n2/d2, n1/d1, and n0/d0 are our successively more accurate
* approximations of the rational. They are, respectively,
* the current, previous, and two prior iterations of it.
*
* a is current term of the continued fraction.
*/
unsigned long n, d, n0, d0, n1, d1, n2, d2;
n = given_numerator;
d = given_denominator;
n0 = d1 = 0;
n1 = d0 = 1;
for (;;) {
unsigned long dp, a;
if (d == 0)
break;
/* Find next term in continued fraction, 'a', via
* Euclidean algorithm.
*/
dp = d;
a = n / d;
d = n % d;
n = dp;
/* Calculate the current rational approximation (aka
* convergent), n2/d2, using the term just found and
* the two prior approximations.
*/
n2 = n0 + a * n1;
d2 = d0 + a * d1;
/* If the current convergent exceeds the maxes, then
* return either the previous convergent or the
* largest semi-convergent, the final term of which is
* found below as 't'.
*/
if ((n2 > max_numerator) || (d2 > max_denominator)) {
unsigned long t = ULONG_MAX;
if (d1)
t = (max_denominator - d0) / d1;
if (n1)
t = min(t, (max_numerator - n0) / n1);
/* This tests if the semi-convergent is closer than the previous
* convergent. If d1 is zero there is no previous convergent as this
* is the 1st iteration, so always choose the semi-convergent.
*/
if (!d1 || 2u * t > a || (2u * t == a && d0 * dp > d1 * d)) {
n1 = n0 + t * n1;
d1 = d0 + t * d1;
}
break;
}
n0 = n1;
n1 = n2;
d0 = d1;
d1 = d2;
}
*best_numerator = n1;
*best_denominator = d1;
}
EXPORT_SYMBOL(rational_best_approximation);