任意の有理数に対する「数当てゲーム」？

質問

以前、面接の質問で次のようなことを言われたことがあります。

私は正の整数nを考えているのですが、それをO(lg n)クエリで推測できるアルゴリズムを考えてください。 各クエリはあなたが選んだ数であり、私は "lower," "higher," または "correct." のいずれかに答えます。

この問題は、2の累乗をnを超えるものが見つかるまでリストアップし、その範囲に対して標準的な二項探索を実行するという、修正二項探索によって解決することができます。 これに関して私がとてもクールだと思うのは、特定の数に対してブルート フォースよりも高速に無限の空間を検索できることです。

しかし、私が抱いている疑問は、この問題のわずかな修正版です。 正の整数を選ぶ代わりに、私が 任意の有理数 を選んだとします。 私の質問は、私が選んだ有理数を最も効率的に決定するために、どのようなアルゴリズムが使えるか、ということです。

今現在、私が持っている最も良い解決策は、暗黙のうちに シュテルン-ブロコット木 というすべての有理数上の二項探索木を暗黙的に歩くことによって p/q を見つけることができます。 しかし、私は整数の場合に得られた実行時間に近いものを得たいと思いました。おそらくO(lg (p + q)) か O(lg pq) のようなものです。 このような実行時間を得るための方法を知っている人はいますか？

私は最初、区間 [0, 1] の標準的なバイナリ検索を使用することを考えましたが、これは非反復バイナリ表現を持つ有理数のみを見つけ、ほとんどすべての有理数を見逃します。 また、有理数を列挙する他の方法を使用することも考えましたが、大きい/等しい/小さいの比較だけで、この空間を検索する方法を見つけることができないようです。

どのように解決するのですか?

さて、以下が私の答えです。 連続分数 を単独で使っての答えです。

まず、ここで用語を整理しておきましょう。

X＝p／qを未知数分数とする。

Q(X,p/q) = sign(X - p/q)を問合せ関数とします。これが0であれば、数を推測したことになり、+/-1であれば、間違いの符号がわかります。

は 連続分数に対する従来の表記法 は A = [a <サブ 0 ; a ₁ , a ₂ , a ₃ , ... a _k ]

= a <サブ 0 + 1/(a ₁ + 1/(a ₂ + 1/(a ₃ + 1/( ... + 1/a <サブ k ) ... )))

0 < p/q < 1の場合、以下のアルゴリズムに従うことにします。

1. Y = 0 = [ 0 ], Z = 1 = [ 1 ], k = 0 を初期化する。

2. 外周ループ : その 前提条件 がそれにあたります。

   • YとZはk+1項の継続分数で、最後の要素以外は同じであり、1だけ異なるので、Y = [y]となる。 ₀ ; y ₁ , y ₂ , y ₃ , ... y <サブ k とし，Z = [y ₀ ; y ₁ , y ₂ , y ₃ , ... y <サブ k + 1]

   • (-1) ^k (Y-X) < 0 < (-1) ^k (Z-X)、またはより簡単に言うと、kが偶数の場合はY < X < Z、kが奇数の場合はZ < X < Yです。

3. 数の値を変えずに、継続分数の次数を1段階拡張する。一般に、最後の項がyの場合 _k であり、y _k + 1 とすれば，これを [... y _k , y _k+1 =∞] となり，[... y _k , z _k+1 =1]. ここでkを1だけ増やす。

4. 内部ループ : これは@templatetypedefさんの整数に関するインタビューの質問と本質的に同じです。二相二値探索を行い、近づける。

5. 内部ループ1 : y <サブ k = ∞, z <サブ k = a であり、X は Y と Z の間にある。

6. Zの最終項を2倍にする。M = Z を計算するが、m <サブ k = 2*a = 2*z <サブ k .

7. 未知数を問い合わせる：q = Q(X,M).

8. q = 0の場合、私たちは答えを持っているので、ステップ17に進みます。

9. qとQ(X,Y)の符号が逆なら、XがYとMの間にあることを意味するので、Z＝Mとし、手順5へ進む。

10. そうでなければY＝Mと設定し、次のステップに進みます。

11. インナーループ2。 y <サブ k = b, z <サブ k = a であり、X は Y と Z の間にある。

12. aとbが1だけ異なる場合、YとZを入れ替え、手順2へ進みます。

13. 二項探索を行う：Mを計算し、ここでmは <サブ k = floor((a+b)/2, query q = Q(X,M))とする。

14. q = 0なら終了、ステップ17に進みます。

15. qとQ(X,Y)の符号が逆なら、XがYとMの間にあることを意味するので、Z＝Mとし、ステップ11に進む。

16. そうでなければ、qとQ(X,Z)が反対の符号を持つ場合、XがZとMの間にあることを意味するので、Y=Mと設定して手順11に進みます。

17. 完了です。X = Mです。

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!

Mを計算する各ステップで、区間の範囲は減少します。各ステップで区間が少なくとも1/sqrt(5)のファクターで減少することを証明するのはおそらくかなり簡単で、それはこのアルゴリズムがO(log q)ステップであることを示すことになるでしょう。

これはtemplatetypedefのオリジナルのインタビューの質問と組み合わされ、次のように適用することができることに注意してください。 任意の 有理数 p/q に対して適用することができることに注意してください。

次の機会があれば、このアルゴリズムを実装したPythonプログラムを投稿します。

編集 : また、各ステップで継続分数を計算する必要がないことにも注意してください (これは O(k) となります。継続分数の部分近似があり、前のステップから次のステップを O(1) で計算することができます)。

2を編集 : 部分近似の再帰的定義。

もしA _k = [a ₀ ; a ₁ , a ₂ , a ₃ , ... a _k ] = p _k /q <サブ k とすると、p _k = a _k p <サブ k-1 + p _k-2 であり，q _k = a _k q <サブ k-1 + q _k-2 . (出典 Niven & Zuckerman, 4th ed, Theorems 7.3-7.5. 参照 ウィキペディア )

例 [0] = 0/1 = p <サブ 0 /q <サブ 0 , [0; 7] = 1/7 = p ₁ /q ₁ ということで、[0; 7, 16] = (16*1+0)/(16*7+1) = 16/113 = p ₂ /q ₂ .

これは、2つの継続分数YとZが最後の項を除いて同じであり、最後の項を除いた継続分数をpとすると _k-1 /q _k-1 とすれば，Y = (y _k p _k-1 + p _k-2 ) / (y _k q <サブ k-1 + q _k-2 ) とし、Z = (z _k p _k-1 + p _k-2 ) / (z _k q <サブ k-1 + q _k-2 ). ここから、このアルゴリズムによって生成される各小間隔で、｜Y-Z｜が少なくとも1/sqrt(5)の要因によって減少することを示すことができるはずですが、代数は今のところ私を超えているようです :-()

これが私の Python プログラムです。

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]

の出力例と 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

Pythonは最初から大整数の計算を扱っており、このプログラムでは（区間計算を除いて）整数計算しか使っていないので、任意の有理数に対して動作するはずです。

3を編集 : O(log^2 q)ではなく、O(log q)であることの証明の概要です。

まず、有理数が見つかるまでのステップ数nは _k は新しい連続分数項ごとに まさに 2b(a_k)-1 ここで b(a_k) は a_k = ceil(log2(a_k)) を表すのに必要なビット数です： バイナリサーチの "net" を広げるには b(a_k) ステップ、狭くするには b(a_k)-1 ステップとなります）。上の例を見てください。ステップ数は常に1、3、7、15などであることに注意してください。

ここで、漸化式関係 q _k = a _k q <サブ k-1 + q _k-2 と帰納法を用いて目的の結果を証明します。

このように述べると、Nの後のqの値は _k = sum(n _k ) k番目の項に到達するのに必要なステップ数が最小となる： q >= A*2 ^cN である（つまり、反転させると、ステップ数Nは <= (1/c) * log ₂ (q/A) = O(log q)です)。

基本的なケース

• k=0：q=1、N=0なのでq >=2 ^N
• k=1：N=2b-1ステップの場合、q=a ₁ >= 2 ^b-1 = 2 ^(N-1)/2 = 2 ^N/2 /sqrt(2)です。

これは、A = 1、c = 1/2が望ましい境界を提供できることを意味します。現実には、q は ではない は各項を2倍にする（反例。[0; 1, 1, 1, 1, 1] は phi = (1+sqrt(5))/2 の成長因子を持つ）ので、c = 1/4 を使いましょう。

帰納法です。

• 項k, qに対して <サブ k = a _k q <サブ k-1 + q _k-2 . ここでも，n _k = 2b-1 のステップを必要とするため、この項では _k >= 2 ^b-1 = 2 ^{(n _k -1)/2} .

  ということは _k q <サブ k-1 >= 2 ^{(N <サブ k -1)/2} * q <サブ k-1 >= 2 ^{(n <サブ k -1)/2} * A*2 ^{N _k-1 /4} = A*2 ^{N _k /4} /sqrt(2)*2 ^{n _k /4} .

アーッ -- ここで大変なのは、もし _k = 1 の場合、q はその一項ではあまり増加しない可能性があり、q を _k-2 を使う必要がありますが、これは q よりもずっと小さいかもしれません。 _k-1 .