algorithm - The "guess the number" game for arbitrary rational numbers?

Question

I once got the following as an interview question:

I'm thinking of a positive integer n. Come up with an algorithm that can guess it in O(lg n) queries. Each query is a number of your choosing, and I will answer either "lower," "higher," or "correct."

This problem can be solved by a modified binary search, in which you listing powers of two until you find one that exceeds n, then run a standard binary search over that range. What I think is so cool about this is that you can search an infinite space for a particular number faster than just brute-force.

The question I have, though, is a slight modification of this problem. Instead of picking a positive integer, suppose that I pick an arbitrary rational number between zero and one. My question is: what algorithm can you use to most efficiently determine which rational number I've picked?

Right now, the best solution I have can find p/q in at most O(q) time by implicitly walking the Stern-Brocot tree, a binary search tree over all the rationals. However, I was hoping to get a runtime closer to the runtime that we got for the integer case, maybe something like O(lg (p + q)) or O(lg pq). Does anyone know of a way to get this sort of runtime?

I initially considered using a standard binary search of the interval [0, 1], but this will only find rational numbers with a non-repeating binary representation, which misses almost all of the rationals. I also thought about using some other way of enumerating the rationals, but I can't seem to find a way to search this space given just greater/equal/less comparisons.

Answer 1

Okay, here's my answer using continued fractions alone.

First let's get some terminology here.

Let X = p/q be the unknown fraction.

Let Q(X,p/q) = sign(X - p/q) be the query function: if it is 0, we've guessed the number, and if it's +/- 1 that tells us the sign of our error.

The conventional notation for continued fractions is A = [a₀; a₁, a₂, a₃, ... a_k]

= a₀ + 1/(a₁ + 1/(a₂ + 1/(a₃ + 1/( ... + 1/a_k) ... )))

---

We'll follow the following algorithm for 0 < p/q < 1.

1. Initialize Y = 0 = [ 0 ], Z = 1 = [ 1 ], k = 0.

2. Outer loop: The preconditions are that:

   • Y and Z are continued fractions of k+1 terms which are identical except in the last element, where they differ by 1, so that Y = [y₀; y₁, y₂, y₃, ... y_k] and Z = [y₀; y₁, y₂, y₃, ... y_k + 1]

   • (-1)^k(Y-X) < 0 < (-1)^k(Z-X), or in simpler terms, for k even, Y < X < Z and for k odd, Z < X < Y.

3. Extend the degree of the continued fraction by 1 step without changing the values of the numbers. In general, if the last terms are y_k and y_k + 1, we change that to [... y_k, y_k+1=∞] and [... y_k, z_k+1=1]. Now increase k by 1.

4. Inner loops: This is essentially the same as @templatetypedef's interview question about the integers. We do a two-phase binary search to get closer:

5. Inner loop 1: y_k = ∞, z_k = a, and X is between Y and Z.

6. Double Z's last term: Compute M = Z but with m_k = 2*a = 2*z_k.

7. Query the unknown number: q = Q(X,M).

8. If q = 0, we have our answer and go to step 17 .

9. If q and Q(X,Y) have opposite signs, it means X is between Y and M, so set Z = M and go to step 5.

10. Otherwise set Y = M and go to the next step:

11. Inner loop 2. y_k = b, z_k = a, and X is between Y and Z.

12. If a and b differ by 1, swap Y and Z, go to step 2.

13. Perform a binary search: compute M where m_k = floor((a+b)/2, and query q = Q(X,M).

14. If q = 0, we're done and go to step 17.

15. If q and Q(X,Y) have opposite signs, it means X is between Y and M, so set Z = M and go to step 11.

16. Otherwise, q and Q(X,Z) have opposite signs, it means X is between Z and M, so set Y = M and go to step 11.

17. Done: X = M.

A concrete example for X = 16/113 = 0.14159292

Y = 0 = [0], Z = 1 = [1], k = 0

k = 1:
Y = 0 = [0; &#8734;] < X, Z = 1 = [0; 1] > X, M = [0; 2] = 1/2 > X.
Y = 0 = [0; &#8734;], Z = 1/2 = [0; 2], M = [0; 4] = 1/4 > X.
Y = 0 = [0; &#8734;], Z = 1/4 = [0; 4], M = [0; 8] = 1/8 < X.
Y = 1/8 = [0; 8], Z = 1/4 = [0; 4], M = [0; 6] = 1/6 > X.
Y = 1/8 = [0; 8], Z = 1/6 = [0; 6], M = [0; 7] = 1/7 > X.
Y = 1/8 = [0; 8], Z = 1/7 = [0; 7] 
  --> the two last terms differ by one, so swap and repeat outer loop.

k = 2:
Y = 1/7 = [0; 7, &#8734;] > X, Z = 1/8 = [0; 7, 1] < X,
    M = [0; 7, 2] = 2/15 < X
Y = 1/7 = [0; 7, &#8734;], Z = 2/15 = [0; 7, 2],
    M = [0; 7, 4] = 4/29 < X
Y = 1/7 = [0; 7, &#8734;], Z = 4/29 = [0; 7, 4], 
    M = [0; 7, 8] = 8/57 < X
Y = 1/7 = [0; 7, &#8734;], Z = 8/57 = [0; 7, 8],
    M = [0; 7, 16] = 16/113 = X 
    --> done!

At each step of computing M, the range of the interval reduces. It is probably fairly easy to prove (though I won't do this) that the interval reduces by a factor of at least 1/sqrt(5) at each step, which would show that this algorithm is O(log q) steps.

Note that this can be combined with templatetypedef's original interview question and apply towards any rational number p/q, not just between 0 and 1, by first computing Q(X,0), then for either positive/negative integers, bounding between two consecutive integers, and then using the above algorithm for the fractional part.

When I have a chance next, I will post a python program that implements this algorithm.

edit: also, note that you don't have to compute the continued fraction each step (which would be O(k), there are partial approximants to continued fractions that can compute the next step from the previous step in O(1).)

edit 2: Recursive definition of partial approximants:

If A_k = [a₀; a₁, a₂, a₃, ... a_k] = p_k/q_k, then p_k = a_kp_k-1 + p_k-2, and q_k = a_kq_k-1 + q_k-2. (Source: Niven & Zuckerman, 4th ed, Theorems 7.3-7.5. See also Wikipedia)

Example: [0] = 0/1 = p₀/q₀, [0; 7] = 1/7 = p₁/q₁; so [0; 7, 16] = (16*1+0)/(16*7+1) = 16/113 = p₂/q₂.

This means that if two continued fractions Y and Z have the same terms except the last one, and the continued fraction excluding the last term is p_k-1/q_k-1, then we can write Y = (y_kp_k-1 + p_k-2) / (y_kq_k-1 + q_k-2) and Z = (z_kp_k-1 + p_k-2) / (z_kq_k-1 + q_k-2). It should be possible to show from this that |Y-Z| decreases by at least a factor of 1/sqrt(5) at each smaller interval produced by this algorithm, but the algebra seems to be beyond me at the moment. :-(

Here's my Python program:

import math

# Return a function that returns Q(p0/q0,p/q) 
#   = sign(p0/q0-p/q) = sign(p0q-q0p)*sign(q0*q)
# If p/q < p0/q0, then Q() = 1; if p/q < p0/q0, then Q() = -1; otherwise Q()=0.
def makeQ(p0,q0):
  def Q(p,q):
    return cmp(q0*p,p0*q)*cmp(q0*q,0)
  return Q

def strsign(s):
  return '<' if s<0 else '>' if s>0 else '=='

def cfnext(p1,q1,p2,q2,a):
  return [a*p1+p2,a*q1+q2]

def ratguess(Q, doprint, kmax):
# p2/q2 = p[k-2]/q[k-2]
  p2 = 1
  q2 = 0
# p1/q1 = p[k-1]/q[k-1]
  p1 = 0
  q1 = 1
  k = 0
  cf = [0]
  done = False
  while not done and (not kmax or k < kmax):
    if doprint:
      print 'p/q='+str(cf)+'='+str(p1)+'/'+str(q1)
# extend continued fraction
    k = k + 1
    [py,qy] = [p1,q1]
    [pz,qz] = cfnext(p1,q1,p2,q2,1)
    ay = None
    az = 1
    sy = Q(py,qy)
    sz = Q(pz,qz)
    while not done:
      if doprint:
        out = str(py)+'/'+str(qy)+' '+strsign(sy)+' X '
        out += strsign(-sz)+' '+str(pz)+'/'+str(qz)
        out += ', interval='+str(abs(1.0*py/qy-1.0*pz/qz))
      if ay:
        if (ay - az == 1):
          [p0,q0,a0] = [pz,qz,az]
          break
        am = (ay+az)/2
      else:
        am = az * 2
      [pm,qm] = cfnext(p1,q1,p2,q2,am)
      sm = Q(pm,qm)
      if doprint:
        out = str(ay)+':'+str(am)+':'+str(az) + '   ' + out + ';  M='+str(pm)+'/'+str(qm)+' '+strsign(sm)+' X '
        print out
      if (sm == 0):
        [p0,q0,a0] = [pm,qm,am]
        done = True
        break
      elif (sm == sy):
        [py,qy,ay,sy] = [pm,qm,am,sm]
      else:
        [pz,qz,az,sz] = [pm,qm,am,sm]     

    [p2,q2] = [p1,q1]
    [p1,q1] = [p0,q0]    
    cf += [a0]

  print 'p/q='+str(cf)+'='+str(p1)+'/'+str(q1)
  return [p1,q1]

and a sample output for ratguess(makeQ(33102,113017), True, 20):

p/q=[0]=0/1
None:2:1   0/1 < X < 1/1, interval=1.0;  M=1/2 > X 
None:4:2   0/1 < X < 1/2, interval=0.5;  M=1/4 < X 
4:3:2   1/4 < X < 1/2, interval=0.25;  M=1/3 > X 
p/q=[0, 3]=1/3
None:2:1   1/3 > X > 1/4, interval=0.0833333333333;  M=2/7 < X 
None:4:2   1/3 > X > 2/7, interval=0.047619047619;  M=4/13 > X 
4:3:2   4/13 > X > 2/7, interval=0.021978021978;  M=3/10 > X 
p/q=[0, 3, 2]=2/7
None:2:1   2/7 < X < 3/10, interval=0.0142857142857;  M=5/17 > X 
None:4:2   2/7 < X < 5/17, interval=0.00840336134454;  M=9/31 < X 
4:3:2   9/31 < X < 5/17, interval=0.00379506641366;  M=7/24 < X 
p/q=[0, 3, 2, 2]=5/17
None:2:1   5/17 > X > 7/24, interval=0.00245098039216;  M=12/41 < X 
None:4:2   5/17 > X > 12/41, interval=0.00143472022956;  M=22/75 > X 
4:3:2   22/75 > X > 12/41, interval=0.000650406504065;  M=17/58 > X 
p/q=[0, 3, 2, 2, 2]=12/41
None:2:1   12/41 < X < 17/58, interval=0.000420521446594;  M=29/99 > X 
None:4:2   12/41 < X < 29/99, interval=0.000246366100025;  M=53/181 < X 
4:3:2   53/181 < X < 29/99, interval=0.000111613371282;  M=41/140 < X 
p/q=[0, 3, 2, 2, 2, 2]=29/99
None:2:1   29/99 > X > 41/140, interval=7.21500721501e-05;  M=70/239 < X 
None:4:2   29/99 > X > 70/239, interval=4.226364059e-05;  M=128/437 > X 
4:3:2   128/437 > X > 70/239, interval=1.91492009996e-05;  M=99/338 > X 
p/q=[0, 3, 2, 2, 2, 2, 2]=70/239
None:2:1   70/239 < X < 99/338, interval=1.23789953207e-05;  M=169/577 > X 
None:4:2   70/239 < X < 169/577, interval=7.2514738621e-06;  M=309/1055 < X 
4:3:2   309/1055 < X < 169/577, interval=3.28550190148e-06;  M=239/816 < X 
p/q=[0, 3, 2, 2, 2, 2, 2, 2]=169/577
None:2:1   169/577 > X > 239/816, interval=2.12389981991e-06;  M=408/1393 < X 
None:4:2   169/577 > X > 408/1393, interval=1.24415093544e-06;  M=746/2547 < X 
None:8:4   169/577 > X > 746/2547, interval=6.80448470014e-07;  M=1422/4855 < X 
None:16:8   169/577 > X > 1422/4855, interval=3.56972657711e-07;  M=2774/9471 > X 
16:12:8   2774/9471 > X > 1422/4855, interval=1.73982239227e-07;  M=2098/7163 > X 
12:10:8   2098/7163 > X > 1422/4855, interval=1.15020646951e-07;  M=1760/6009 > X 
10:9:8   1760/6009 > X > 1422/4855, interval=6.85549088053e-08;  M=1591/5432 < X 
p/q=[0, 3, 2, 2, 2, 2, 2, 2, 9]=1591/5432
None:2:1   1591/5432 < X < 1760/6009, interval=3.06364213998e-08;  M=3351/11441 < X 
p/q=[0, 3, 2, 2, 2, 2, 2, 2, 9, 1]=1760/6009
None:2:1   1760/6009 > X > 3351/11441, interval=1.45456726663e-08;  M=5111/17450 < X 
None:4:2   1760/6009 > X > 5111/17450, interval=9.53679318849e-09;  M=8631/29468 < X 
None:8:4   1760/6009 > X > 8631/29468, interval=5.6473816179e-09;  M=15671/53504 < X 
None:16:8   1760/6009 > X > 15671/53504, interval=3.11036635336e-09;  M=29751/101576 > X 
16:12:8   29751/101576 > X > 15671/53504, interval=1.47201634215e-09;  M=22711/77540 > X 
12:10:8   22711/77540 > X > 15671/53504, interval=9.64157420569e-10;  M=19191/65522 > X 
10:9:8   19191/65522 > X > 15671/53504, interval=5.70501257346e-10;  M=17431/59513 > X 
p/q=[0, 3, 2, 2, 2, 2, 2, 2, 9, 1, 8]=15671/53504
None:2:1   15671/53504 < X < 17431/59513, interval=3.14052228667e-10;  M=33102/113017 == X

Since Python handles biginteger math from the start, and this program uses only integer math (except for the interval calculations), it should work for arbitrary rationals.