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

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 2

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 3

There are two ways to do that. One is to stick with some form of tree (graph) traversal algorithm -- then what you need is to write an algorithm which traverses Stern-Brocot tree. It is a binary tree but you will need to check at every step if the fraction that represents your current node is in $F_M$ (Farey series of order M, which is what you defined as $R(M)$) -- Stern-Brocot tree of order M will have lots of nodes/fractions with denominators greater than M. The way to write such traversal code is to use continued fraction representations as node data and then proceed using the fact that continued fraction $[a_0; a_1, \ldots, a_k]$ has two children $[a_0; a_1, \ldots, a_k+1]$ and $[a_0; a_1, \ldots, a_k-1, 2]$, etc. Start with the root $1/2 = [0;2]$, slide left-down to $1/M$ etc -- I assume you don't want to use recursion.

The other way is much simpler but less exciting :), and it might not be available to you if you only have one member of $F_M$ (however if you are simply looking for a way to quickly write out Farey series then you are fine). It uses not one but two last numbers -- current one $p/q$ and previous one $u/v$. Then we know that the next number $s/t$ produces $p/q$ as its mediant with $u/v$, ie $$ \frac{p}{q} = \frac{u+s}{v+t}\quad\Rightarrow\quad u+s = kp, v+t = kq $$ for some natural number $k$. This gives us $s = kp-u, t = kq-v$, we just need to find $k$. Obviously the larger $k$ the closer $s/t$ is to $p/q$, so we need the maximum possible value of $k$. As an exercise prove that $k = \lfloor(M+v)/q\rfloor$ is the answer (sorry, educator's habit)...

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

Answer 6

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?