// 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> | |

#include <linux/module.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); | |

MODULE_LICENSE("GPL v2"); |