Internet Problem Solving Contest

IPSC 2005

Solution to Problem E – Exact? Approximate!

First of all, we hope that no calculator was (or will be) harmed as a consequence of this contest.

One possible way to solve the easy input was to increase the denominator until we find an integer numerator such that the resulting fraction is close enough to the input value. But solving the hard input this way required a stronger mathematical background.

Let X/Y be our input number. (X are the digits after the decimal point, Y=10D.) How can we tell whether a fraction A/B approximates it well? It's simple. The value A/B rounded to D decimal places should give X/Y. Thus the exact value has to be from the interval < X/Y - 0.5*(1/Y) , X/Y + 0.5*(1/Y) ) – or equivalently, from the interval < (10X-5) / 10Y , (10X+5) / 10Y ).

Now we need to know how to compare two fractions. This can be done using integer arithmetics only – note that for positive numbers A, B, C, D we have A/B < C/D iff A*D < B*C .

Now consider a structure called the Stern-Brocot tree . Truncated at some level and recursively written in left-center-right manner, it guarantees that it lists rationals in increasing order, with the denominator bounded by some constant. (Such representation of Stern-Brocot tree is called the Farey sequence).

Thus, if we search the tree from its root for a fraction approximating our input number X/Y, we can be sure that the first fraction that approximates X/Y is the one we are looking for (because all the following fractions will have larger denominators).

Looking at the input, we realize that ordinary integer or floating-point types in C or Pascal don't have sufficient precision. Therefore, we either implement our own arithmetics (integer addition, multiplication and comparison should be enough), or use some tool with arbitrary precision numbers, such as the linux calculator bc.

We can now code our first attempt. We will do something like a binary search – we will keep the lower bound A/B and the upper bound C/D of some interval in which we currently search. The initial bounds will be 0/1 and 1/1.

In each step, we compute the "middle" fraction E/F as (A+C)/(B+D). Note that this will always be the fraction with the smallest denominator in the current interval. (This is a property of the Stern-Brocot tree mentioned above.)

If E/F approximates X/Y, we are finished. Otherwise, we will compare these two fractions and change the lower or the upper bound to E/F.

Unfortunately, the program, although correct, isn't fast enough to compute larger inputs. Therefore, we must optimize it. But how, where can it be optimized?

If we watch the path we take in the tree for a given input, we can realize that sometimes we take many successive steps to one side. Let us consider the input 0.94723577432. We will first take 17 steps to the right before taking the first step to the left. If we were able to take all these steps at once, we would save many comparisons and multiplications of large integers. This is definitely worth trying.

Therefore, when computing the middle fraction, we try to do as many steps at once as possible. For example, when we go to the right, only the lower bound changes. Therefore we know, that after K steps, we will still have the same upper bound C/D. The middle fraction after K steps to the right would be (A+K*C)/(B+K*D). We now have to find the largest K such that this fraction is still less than X/Y (otherwise we would be forced to take a left step sooner). By solving the inequality (A+K*C)/(B+K*D) < X/Y we get that the largest integer K satisfying it is K=ceiling( (B*X - A*Y) / (C*Y - D*X) ) - 1.

Note that if K=0, let K=1 – we really should make at least one step, even if we exceed X/Y. The other case (going to the left) is handled symetrically, read our code for clarification.

Now our program simply takes the value E/F = (A+K*C)/(B+K*D) as the new "midpoint" and we are done. Wait... Is that really so? Not exactly. Consider the case when E/F approximates X/Y well. Now it doesn't have to be true that the denominator is as small as possible – what if one of the fractions we skipped was already close enough?

Luckily, if both (A+K*C)/(B+K*D) and (A+L*C)/(B+L*D) approximate X/Y well, so do all the fractions (A+q*C)/(B+q*D) with q between K and L. Thus we may use binary search to find the smallest good value of L. (There were some inputs where a linear search would take ages to finish.)

It can be shown that after a finite number of steps (8 should be a convenient upper bound) the denominators of both ends of the current interval will at least double. Also, X/Y is a good approximation of X/Y. Thus the time complexity of finding the first solution is polynomial in the size of the input. Thanks to the binary search, so is finding the smallest solution from the first found one.

Note: There is a function "Rationalize" in Mathematica, which converts rational numbers to fractions. First of all note that using Mathematica was forbidden. Second of all, the function in Mathematica sometimes gives a slightly non-optimal result. Thus even if you did use it, your answer wouldn't be accepted.