How to compute next/previous representable rational number?

Question

An (approximate) non-negative rational number representation is a pair of natural numbers each not greater than some fixed limit M (and of course denominator being non-zero).

With this condition there is finite number of representable rational numbers. This means that for each such number we can name previous and next number in the set (of course except for smallest 0 and largest 1). How to compute them?

In other words lets have a set of numbers

$$ R(M)=\{\frac{p}{q} : p,q\in\mathbb{N},q\neq0,p\leq q\leq M\} $$

where $M\in\mathbb{N}_+$

and given number $\frac{p_1}{q_1}\in R(M)$.

I want to find numbers $\frac{p_S}{q_S}\in R(M)$ and $\frac{p_L}{q_L}\in R(M)$ (if they exist) such that

$$ \frac{p_S}{q_S}<\frac{p_1}{q_1} \wedge \neg\exists_{\frac{p}{q}\in R(M)}\frac{p_S}{q_S}<\frac{p}{q}<\frac{p_1}{q_1} $$

and smilarly

$$ \frac{p_L}{q_L}>\frac{p_1}{q_1} \wedge \neg\exists_{\frac{p}{q}\in R(M)}\frac{p_1}{q_1}<\frac{p}{q}<\frac{p_L}{q_L} $$

Answer 1

All this needs is a single modular inversion, with complexity $O(\log M \log \log M)$ (or something like that).

We are given $p/q$ in lowest terms, with $p, q \in \mathbb{N}$ and $q > 0$. We want the nearest neighbours of $p/q$, $a/b < p/q < c/d$, subject to $a,b,c,d \in \mathbb{N}$ and $c,d > 0$. Assume for the moment that $p < q$.

Consider $c/d$. It is the element after $p/q$ in the Farey sequence $F_n$, which means that $cq - dp = 1$. In particular, $dp+1 = 0 \pmod q$, i.e. $d = -1/p \pmod q$. So let

$r = 1/p \pmod q$

Now make $d$ as big as possible, subject to (i) $d = -r \pmod q$ and (ii) $d \leq n$. Then just put $c = (dp+1)/q$.

Similarly, to calculate $a/b$, make $b$ as big as possible subject to (i) $b = +r \pmod q$ and (ii) $b \leq n$. Then put $a = (bp-1)/q$.

An example Suppose $p/q = 28/61$ and $n=100$. Then the inverse of $28 \pmod{61}$ is $r=24$. So $d = 37 \pmod{61}$, and we can take $d=98$, with $c=(98.28+1)/61 = 45$. And $b = 24 \pmod{61}$, so we can take $b = 85$, with $a = (85.28-1)/61=39$. So we get:

$39/85 < 28/61 < 45/98$

We only have to modify this slightly to handle the case $p > q$: instead of computing the largest $d \leq n$ that satisfies the modulo requirement, we find the corresponding modulo requirement for $c$, and then look for the largest such $c \leq n$. You can fill in the details for yourself.

Answer 2

Check out the next link: http://en.wikipedia.org/wiki/Farey_sequence . In general, Farey sequence is dealing exactly with that type of question.

Answer 3

So far I have found algorithm iterating over denominators (inspired by mjqxxxx's comment). Thus it is of $O(M)$ time complexity (however each iteration is short and simple).

In Farey sequence for two consecutive fractions $\frac{p_1}{q_1}<\frac{p_L}{q_L}$ it holds that $p_L q_1 - p_1 q_L=1$. Thus $p_L=\frac{p_1 q_L + 1}{q_1}$. We have $\frac{p_1}{q_1}$ and $M$ and we want to know $\frac{p_L}{q_L}$.

Since $1\leq q_L\leq M$ we can now simply iterate over all natural numbers in range $[1;M]$ and for each see whether $\frac{p_1 q_L + 1}{q_1}$ is a natural number. If so than we have the desired fraction $\frac{p_L}{q_L}$.

It may happen that this will be solved for more than one value in that range. This means that depending on $M$ the solution are different fractions. Thus of interest is the last (I think - I haven't prove that) value and it would be better to iterate from $M$ to $1$ rather than the other way.

(Finding previous fraction can be done in similar way using $p_S=\frac{p_1 q_S – 1}{q_1}$.)

It remains an open question whether indeed the largest $q_L$/$q_S$ which solves the equation in natural numbers is the solution. Also this is still $O(M)$ so maybe something faster can be found?

Answer 4

As I mentioned in a closely related question, this has a very beautiful geometrical interpretation employing Pick's Area Theorem: $\rm\: a/b < c/d\:$ are two adjacent terms in a Farey sequence iff the triangle $\rm T$ with vertices $\rm (0,0), u = (b,a), v = (d,c)$ contains no lattice points except for vertices. By Pick's formula this is true iff $\rm\ area\ T = $ #interior_points $\ +\ 1/2\ $ #boundary_points $\:-\ 1\ =$ $\ 0 + 3/2 - 1 = 1/2\:.$ But by basic analytic geometry $\rm\: area\ T\: =\: |det(u,v)|/2 = |ad-bc|/2\:.\:$ Hence $\rm\:a/b,\:c/d\:$ are meighbors in some Farey sequence iff $\rm\:|ad-bc| = 1\:.\:$

Therefore the problem reduces to solving the linear Diophantine equation $\rm\:ax+by = 1\:.\:$ As is well-known, such equations may be solved by employing the extended Euclidean algorithm, which is conveniently implemented by performing the Euclidean algorithm in parallel on an identity-augmented system of equations. For an example worked-out in detail see my post here.

In fact it deserves to be much better known that Pick originally applied his area formula in a similar way to give a beautiful geometric proof of the Bezout linear representation of the gcd.

Answer 5

For the copy-pasters, below I modified the Extended_Euclidean_algorithm from wikibooks to ensure that the gcd is non-negative, and used TonyK's answer to find the Farey neighbours:

def xgcd(a, b):
    """ Extended Euclidean Algorithm: computes (g, x, y),
        such that a*x + b*y = g = gcd(a, b) >= 0. """
    x0, x1, y0, y1 = 0, 1, 1, 0
    while a != 0:
        (q, a), b = divmod(b, a), a
        y0, y1 = y1, y0 - q * y1
        x0, x1 = x1, x0 - q * x1
    return (b, x0, y0) if b >= 0 else (-b, -x0, -y0)

def Farey_neighbours(n, p, q):
    """ Computes the neighbours of p/q in the Farey sequence of order n. """
    assert q != 0
    if q < 0:
        p, q = -p, -q 
    g, r, _ = xgcd(p, q)
    p, q = p // g, q // g
    b = ((n - r) // q) * q + r
    a = (b * p - 1) // q
    d = ((n + r) // q) * q - r
    c = (d * p + 1) // q
    return (a, b), (c, d)

The functions not only reproduce TonyK's example, but also work in case p > q, p < 0, q < 0, or gcd(p,q) > 1:

print(Farey_neighbours(100, 28, 61))        # returns ((39, 85), (45, 98))
print(Farey_neighbours(100, 28+61, 61))     # returns ((39+85, 85), (45+98, 98))
print(Farey_neighbours(100, 28-61, 61))     # returns ((39-85, 85), (45-98, 98))
print(Farey_neighbours(100, 61-28, -61))    # returns ((39-85, 85), (45-98, 98))
print(Farey_neighbours(100, 28*2, 61*2))    # returns ((39, 85), (45, 98))