What is the complexity of a bracketed search using mediants?

Question

I'm trying to estimate the complexity of an algorithm I've written for the Reko decompiler, where I'm trying to "undo" the tranformation done by a compiler to an integer division by a constant $x / n$. The compiler has converted the division into an integer multiplication and a shift: $(x * \lfloor 2^\beta / n \rfloor) >> \beta$, where $\beta$ is the number of bits of the computer's machine word. The resulting constant multiplication is a lot faster than a division in most contemporary architectures, but no longer resembles the original code.

To illustrate: the C statement

y = x / 10;

will be compiled by the Microsoft Visual C++ compiler to following assembly language

mov edx,1999999Ah  ; load 1/10 * 2^32 
imul eax           ; edx:eax = dividend / 10 * 2 ^32 
mov eax,edx        ; eax = dividend / 10

The net result is that the register eax will now have the expected value of y from the source code.

A naive decompiler will decompile the above to

eax = ((long)eax * 0x1999999A) >> 32;

but Reko aims to make the resulting output more legible than that by recovering the constant that was used in the original division.

The algorithm hinted at above is based on the description on this article in Wikipedia. First, the algorithm treats the constant multiplier as the scaled reciprocal $2^\beta / n$. It converts that to a floating point number $2^\beta r_f$ and then scales it down by $2^\beta$ to $r_f$, where $0.0 < r_f <1.0$. The final, expensive step is to bracket the floating point value $r_f$ between two rational numbers $a/b$, $c/d$ (starting with 0/1 and 1/1) and repeatedly compute the mediant $(a + c)/(b + d)$ until some convergence criterion is reached. The result should be the "best" rational approximation $r$ to the reciprocal $r_f$.

Now, if the bracketing was being done with a typical binary search starting between the rationals $0/2^\beta$ and $2^\beta/2^\beta$, and computing the midpoint $(a/b + c/d)/2$, I expect the algorithm to converge in $O(\beta)$ steps. But what is the complexity of the algorithm if the mediant is used instead?

Answer 1

The relation between the Stern–Brocot tree and Farey sequences shows that if $0 < p/q < 1$ and $(p,q) = 1$ (i.e., $p/q$ is a reduced fraction) then $p/q$ is at the $q$th level of the tree. Since the running term of your algorithm is linear at the level in which you terminate, your algorithm takes time $O(q)$, where $p/q$ is the answer; but this is not so helpful.

You haven't specified what your stopping criterion is, but presumably you have some error threshold $\epsilon$. The question therefore is for which Farey sequence is it the case that adjacent terms are at distance at most $2\epsilon$ (and so every point is at distance $\epsilon$ from some point). Using the fact that the distance between adjacent fractions $p_1/q_1,p_2/q_2$ in a Farey sequence is $1/(q_1q_2)$, it is not hard to show that the maximal distance at the $q$th Farey sequence is $1/q$. Therefore if you are aiming at a distance of $\epsilon$, then your algorithm will run in time $O(1/\epsilon)$ in the worst case.

However, "most" adjacent fractions at the $q$th Farey sequence are at distance $O(1/q^2)$, and so on average you would probably expect a running time of $O(\sqrt{1/\epsilon})$ (unfortunately, this average is with respect to the input rather than the algorithm).