任意有理数的“猜数”游戏?

我曾经在面试中问过这样一个问题:

我在考虑一个正整数 n,提出一个算法,可以在 O (lgn)查询中猜测它。每个查询都是您选择的数字,我将回答“低”、“高”或“正确”

这个问题可以通过修改二进制搜索来解决,在这种搜索中,您列出两个幂,直到找到一个超过 n 的幂,然后在这个范围内运行标准的二进制搜索。我觉得这个很酷的地方在于你可以在一个无限的空间里寻找一个特定的数字,而不仅仅是靠蛮力。

不过,我的问题是对这个问题稍作修改。假设我在0和1之间选择一个 arbitrary rational number,而不是选择一个正整数。我的问题是: 你能用什么算法最有效地确定我选择了哪个有理数?

现在,我所知道的最好的解决方案可以通过隐式遍历 Stern-Brocot 树(一个遍历所有有理数的二叉查找树)在最多 o (q)的时间内找到 p/q。然而,我希望得到一个更接近于整数情况下的运行时,比如 O (lg (p + q))或 O (lg pq)。有人知道如何获得这种运行时吗?

我最初考虑使用区间[0,1]的标准二进制搜索,但是这只能找到具有非重复二进制表示的有理数,而且几乎没有找到所有的有理数。我也想过用其他方法来列举这些有理数,但是我似乎找不到一种方法来搜索这个空间,只是给出了更大的/相等的/更小的比较。

9369 次浏览

You can sort rational numbers in a given interval by for example the pair (denominator, numerator). Then to play the game you can

  1. 使用倍增步骤方法查找区间 [0, N]
  2. 给定一个区间 [a, b],求区间中最接近区间中心的分母最小的有理数

然而,这可能仍然是 O(log(num/den) + den)(不确定,这里的早晨太早,使我认为清楚;)

让我们取有理数,以简化的形式,并写出它们的顺序,首先是分母,然后是分子。

1/2, 1/3, 2/3, 1/4, 3/4, 1/5, 2/5, 3/5, 4/5, 1/6, 5/6, ...

我们的第一个猜测是 1/2。然后我们沿着名单继续,直到我们有3个在我们的范围内。然后,我们将采取2个猜测来搜索该名单。然后我们沿着列表继续,直到剩下的范围内有7个。然后我们将采取3个猜测搜索该名单。诸如此类。

n步骤中,我们将介绍第一种 2O(n)可能性,这种可能性正是你所寻找的效率数量级。

Update: People didn't get the reasoning behind this. The reasoning is simple. We know how to walk a binary tree efficiently. There are O(n2) fractions with maximum denominator n. We could therefore search up to any particular denominator size in O(2*log(n)) = O(log(n)) steps. The problem is that we have an infinite number of possible rationals to search. So we can't just line them all up, order them, and start searching.

Therefore my idea was to line up a few, search, line up more, search, and so on. Each time we line up more we line up about double what we did last time. So we need one more guess than we did last time. Therefore our first pass uses 1 guess to traverse 1 possible rational. Our second uses 2 guesses to traverse 3 possible rationals. Our third uses 3 guesses to traverse 7 possible rationals. And our k'th uses k guesses to traverse 2k-1 possible rationals. For any particular rational m/n, eventually it will wind up putting that rational on a fairly big list that it knows how to do a binary search on efficiently.

如果我们做二进制搜索,然后忽略我们所学到的一切,当我们抓取更多的有理数,然后我们会把所有的有理数,包括 m/nO(log(n))的传递。(这是因为到那时,我们将得到一个包含足够多的有理数的传球,这些有理数可以包含到 m/n的所有有理数。)但是每次通过需要更多的猜测,所以这将是 O(log(n)2)的猜测。

然而,实际上我们做的比这要好得多。通过我们的第一次猜测,我们排除了列表中一半的理性因素太大或太小。我们接下来的两个猜测并没有完全把空间分割成四分之一,但是他们也不会离得太远。我们接下来的3次猜测也没有把空间缩小到八分之一,但是也没有太大的差距。诸如此类。当您把它们放在一起时,我确信结果是在 O(log(n))步骤中找到 m/n。虽然我没有证据。

尝试一下: 下面是生成猜测的代码,这样您就可以玩游戏,看看它的效率如何。

#! /usr/bin/python


from fractions import Fraction
import heapq
import readline
import sys


def generate_next_guesses (low, high, limit):
upcoming = [(low.denominator + high.denominator,
low.numerator + high.numerator,
low.denominator, low.numerator,
high.denominator, high.numerator)]
guesses = []
while len(guesses) < limit:
(mid_d, mid_n, low_d, low_n, high_d, high_n) = upcoming[0]
guesses.append(Fraction(mid_n, mid_d))
heapq.heappushpop(upcoming, (low_d + mid_d, low_n + mid_n,
low_d, low_n, mid_d, mid_n))
heapq.heappush(upcoming, (mid_d + high_d, mid_n + high_n,
mid_d, mid_n, high_d, high_n))
guesses.sort()
return guesses


def ask (num):
while True:
print "Next guess: {0} ({1})".format(num, float(num))
if 1 < len(sys.argv):
wanted = Fraction(sys.argv[1])
if wanted < num:
print "too high"
return 1
elif num < wanted:
print "too low"
return -1
else:
print "correct"
return 0


answer = raw_input("Is this (h)igh, (l)ow, or (c)orrect? ")
if answer == "h":
return 1
elif answer == "l":
return -1
elif answer == "c":
return 0
else:
print "Not understood.  Please say one of (l, c, h)"


guess_size_bound = 2
low = Fraction(0)
high = Fraction(1)
guesses = [Fraction(1,2)]
required_guesses = 0
answer = -1
while 0 != answer:
if 0 == len(guesses):
guess_size_bound *= 2
guesses = generate_next_guesses(low, high, guess_size_bound - 1)
#print (low, high, guesses)
guess = guesses[len(guesses)/2]
answer = ask(guess)
required_guesses += 1
if 0 == answer:
print "Thanks for playing!"
print "I needed %d guesses" % required_guesses
elif 1 == answer:
high = guess
guesses[len(guesses)/2:] = []
else:
low = guess
guesses[0:len(guesses)/2 + 1] = []

作为一个尝试的例子,我尝试了101/1024(0.0986328125) ,发现需要20次猜测才能找到答案。我试了0.98765,猜了45次。我尝试了0.0123456789,它需要66次猜测和大约一秒钟来生成它们。(注意,如果你用一个有理数作为参数调用这个程序,它会为你填充所有的猜测。这是一个非常有用的便利。)

我想我找到了一个 O (log ^ 2(p + q))算法。

To avoid confusion in the next paragraph, a "query" refers to when the guesser gives the challenger a guess, and the challenger responds "bigger" or "smaller". This allows me to reserve the word "guess" for something else, a guess for p + q that is not asked directly to the challenger.

这个想法是首先找到 p + q,使用你在问题中描述的算法: 猜一个值 k,如果 k 太小,将其加倍,然后再试一次。然后,一旦你有一个上界和下界,做一个标准的二进制搜索。这需要 O (log (p + q) T)查询,其中 T 是检查猜测所需查询数量的上限。我们去找 T。

我们想用 r + s < = k 检查所有分数 r/s,然后加倍 k 直到 k 是足够大。注意,有一些 o (k ^ 2)分数需要检查给定的 k 值。构建一个包含所有这些值的平衡二叉查找树,然后搜索它来确定 p/q 是否在树中。它使用 O (log k ^ 2) = O (log k)查询来确认 p/q 不在树中。

We will never guess a value of k greater than 2(p + q). Hence we can take T = O(log(p+q)).

When we guess the correct value for k (i.e., k = p + q), we will submit the query p/q to the challenger in the course of checking our guess for k, and win the game.

查询总数为 O (log ^ 2(p + q))。

请记住,(0,1)中的任何有理数都可以表示为不同(正或负)单位分数的有限和。例如,2/3 = 1/2 + 1/6和2/5 = 1/2-1/10。您可以使用它来执行直接的二进制搜索。

我知道了! 你需要做的是使用二分法和 连分数并行搜索。

Bisection will give you a limit toward a specific real number, as represented as a power of two, and continued fractions will take the real number and find the nearest rational number.

如何并行地运行它们如下所示。

在每个步骤中,lu是二分法的上下界。这个想法是,你可以选择将二分法的范围减半,或者添加一个额外的术语作为连分数表示。当 lu的下一个术语与连分数相同时,你就可以进行下一步的连分数搜索,并使用连分数进行查询。否则,使用二分法将范围减半。

Since both methods increase the denominator by at least a constant factor (bisection goes by factors of 2, continued fractions go by at least a factor of phi = (1+sqrt(5))/2), this means your search should be O(log(q)). (There may be repeated continued fraction calculations, so it may end up as O(log(q)^2).)

我们的连分数搜索需要四舍五入到最接近的整数,而不是使用地板(这在下面更清楚)。

上面的例子有点手势,让我们举一个 r = 1/31的具体例子:

  1. L = 0,u = 1,query = 1/2.0不能表示为连分数,所以我们使用二进制搜索直到 l!= 0.

  2. l = 0, u = 1/2, query = 1/4.

  3. L = 0,u = 1/4,query = 1/8.

  4. l = 0, u = 1/8, query = 1/16.

  5. L = 0,u = 1/16,query = 1/32.

  6. L = 1/32,u = 1/16.现在1/l = 32,1/u = 16,它们有不同的 crac 代表,所以保持 bisecting.query = 3/64。

  7. L = 1/32,u = 3/64,query = 5/128 = 1/25.6

  8. L = 1/32,u = 5/128,query = 9/256 = 1/28.4444...

  9. L = 1/32,u = 9/256,query = 17/512 = 1/30.1176... (四舍五入到1/30)

  10. L = 1/32,u = 17/512,query = 33/1024 = 1/31.0303... (四舍五入到1/31)

  11. L = 33/1024,u = 17/512,query = 67/2048 = 1/30.5672... (四舍五入到1/31)

  12. L = 33/1024,u = 67/2048.在这一点上,l 和 u 都有相同的连分数项31,所以现在我们使用连分数猜测。 Query = 1/31.

  13. 查询 = 1/31.

成功了!

对于另一个例子,我们使用16/113(= 355/113-3,其中355/113非常接近 π)。

[待续,我必须去某个地方]


在进一步的思考中,连分数是一种可行的方法,除了确定下一项之外,不要介意二分法。等我回来再说。

这里还有另一种方法来做到这一点。如果有足够的兴趣,我会尽量填写今晚的细节,但我现在不能,因为我有家庭责任。下面是一个应该解释算法的实现的片段:

low = 0
high = 1
bound = 2
answer = -1
while 0 != answer:
mid = best_continued_fraction((low + high)/2, bound)
while mid == low or mid == high:
bound += bound
mid = best_continued_fraction((low + high)/2, bound)
answer = ask(mid)
if -1 == answer:
low = mid
elif 1 == answer:
high = mid
else:
print_success_message(mid)

这就是解释。best_continued_fraction(x, bound)应该做的是找到最后一个近似于 x的连分数,分母最多为 bound。这个算法将采取多对数步骤来完成,并找到非常好的(尽管不总是最好的)近似值。因此,对于每个 bound,我们将得到一个接近于二进制搜索的结果,通过所有可能的分数。有时候,我们不会找到一个特定的分数,直到我们将界增加到比我们应该增加的更远的地方,但是我们不会太远。

这就是多重对数问题的对数数量。

更新: 和完整的工作代码。

#! /usr/bin/python


from fractions import Fraction
import readline
import sys


operations = [0]


def calculate_continued_fraction(terms):
i = len(terms) - 1
result = Fraction(terms[i])
while 0 < i:
i -= 1
operations[0] += 1
result = terms[i] + 1/result
return result


def best_continued_fraction (x, bound):
error = x - int(x)
terms = [int(x)]
last_estimate = estimate = Fraction(0)
while 0 != error and estimate.numerator < bound:
operations[0] += 1
error = 1/error
term = int(error)
terms.append(term)
error -= term
last_estimate = estimate
estimate = calculate_continued_fraction(terms)
if estimate.numerator < bound:
return estimate
else:
return last_estimate


def ask (num):
while True:
print "Next guess: {0} ({1})".format(num, float(num))
if 1 < len(sys.argv):
wanted = Fraction(sys.argv[1])
if wanted < num:
print "too high"
return 1
elif num < wanted:
print "too low"
return -1
else:
print "correct"
return 0


answer = raw_input("Is this (h)igh, (l)ow, or (c)orrect? ")
if answer == "h":
return 1
elif answer == "l":
return -1
elif answer == "c":
return 0
else:
print "Not understood.  Please say one of (l, c, h)"


ow = Fraction(0)
high = Fraction(1)
bound = 2
answer = -1
guesses = 0
while 0 != answer:
mid = best_continued_fraction((low + high)/2, bound)
guesses += 1
while mid == low or mid == high:
bound += bound
mid = best_continued_fraction((low + high)/2, bound)
answer = ask(mid)
if -1 == answer:
low = mid
elif 1 == answer:
high = mid
else:
print "Thanks for playing!"
print "I needed %d guesses and %d operations" % (guesses, operations[0])

与前一个解决方案相比,它在猜测方面似乎更有效,并且操作要少得多。对于101/1024,需要19次猜测和251次操作。对于.98765,它需要27次猜测和623次运算。对于0.0123456789,需要66次猜测和889次操作。对于傻笑和咧嘴一笑,0.0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456787890123456787878901234567778989(这是前一个的10个副本) ,它需要665次猜测和2。

好了,我想我找到了一个 O (lg2 q)算法来解决这个问题,这个算法是基于 Jason S 关于使用连分数的最优秀的见解。我觉得我应该在这里完善算法这样我们就有了一个完整的解决方案,还有一个运行时分析。

算法背后的直觉是,范围内的任何有理数 p/q 都可以写成

A0 + 1/(a1 + 1/(a2 + 1/(a3 + 1/...))

对于 a的适当选择。这叫做 连分数。更重要的是,虽然这些 aabc0可以通过运行分子和分母上的辗转相除法得到。例如,假设我们想以这种方式表示11/14。我们首先注意到14乘以110,所以11/14的粗略近似值是

0 = 0

Now, suppose that we take the reciprocal of this fraction to get 14/11 = 1 3/11. So if we write

0 + (1/1) = 1

我们得到了更接近11/14的结果。现在我们剩下3/11,我们可以再次求倒数,得到11/3 = 3 2/3,所以我们可以考虑

0 + (1/(1 + 1/3)) = 3/4

这是11/14的另一个很好的近似值。现在,我们有2/3,所以考虑倒数,也就是3/2 = 1 1/2。如果我们写

0 + (1/(1 + 1/(3 + 1/)) = 5/6

我们得到了另一个接近11/14的近似值。最后,我们剩下1/2,它的倒数是2/1。如果我们最终写出

0 + (1 / (1 + 1/(3 + 1/(1 + 1/2)))) = (1 / (1 + 1/(3 + 1/(3/2)))) = (1 / (1 + 1/(3 + 2/3)))) = (1 / (1 + 1/(11/3)))) = (1 / (1 + 3/11)) = 1 / (14/11) = 11/14

which is exactly the fraction we wanted. Moreover, look at the sequence of coefficients we ended up using. If you run the extended Euclidean algorithm on 11 and 14, you get that

11 = 0 x 14 + 11—— > a0 = 0 14 = 1 x 11 + 3—— > a1 = 1 11 = 3 x 3 + 2—— > a2 = 3 3 = 2 x 1 + 1—— > a3 = 2

事实证明(使用的数学比我现在知道的还要多!)这不是巧合 p/q 连分数的系数总是由扩展欧几里得算法构成的。这很棒,因为它告诉我们两件事:

  1. 最多可以有 O (lg (p + q))系数,因为辗转相除法总是以这么多步结束
  2. 每个系数最大为{ p,q }。

Given these two facts, we can come up with an algorithm to recover any rational number p/q, not just those between 0 and 1, by applying the general algorithm for guessing arbitrary integers n one at a time to recover all of the coefficients in the continued fraction for p/q. For now, though, we'll just worry about numbers in the range (0, 1], since the logic for handling arbitrary rational numbers can be done easily given this as a subroutine.

作为第一步,假设我们希望找到 a1的最佳值,这样1/a1就尽可能接近 p/q,a1就是一个整数。要做到这一点,我们只需运行我们的算法猜测任意整数,每次取倒数。这样做之后,会发生两件事中的一件。首先,我们可能巧合地发现,对于某个整数 k,p/q = 1/k,在这种情况下,我们完成了。如果不是,我们会发现 p/q 夹在1/(a1-1)和1/a0之间。当我们这样做的时候,我们就可以通过找到 a2使得 p/q 介于1/(a1 + 1/a2)和1/(a1 + 1/(a2 + 1))之间,从而在更深的连分数上开始工作。如果我们能神奇地找到 p/q,那就太好了!否则,我们就在连分数里再往下一层。最终,我们会通过这种方式找到号码,不会花太长时间。每个二进制搜索寻找一个系数最多需要 O (lg (p + q))时间,最多有 O (lg (p + q))级别的搜索,所以我们只需要 O (lg11(p + q))算术运算和探针来恢复 p/q。

我想指出的一个细节是,在搜索时,我们需要跟踪我们是在奇数水平还是在偶数水平,因为当我们在两个连续分数之间夹住 p/q 时,我们需要知道我们所寻找的系数是上部分数还是下部分数。我将在没有证据的情况下声明,对于 i 单数的 a,你需要使用两个数字的上半部分,而对于 a,你甚至需要使用两个数字的下半部分。

I am almost 100% confident that this algorithm works. I'm going to try to write up a more formal proof of this in which I fill in all of the gaps in this reasoning, and when I do I'll post a link here.

感谢大家贡献了使这个解决方案工作所需的见解,特别是 Jason S 建议对连续分数进行二进制搜索。

好的,这是我单独使用 连分数的答案。

首先,让我们在这里得到一些术语。

设 X = p/q 为未知分数。

让 Q (X,p/q) = sign (X-p/q)作为查询函数: 如果是0,我们已经猜出了数字,如果是 +/-1,则告诉我们错误的符号。

连分数的常规记数法是 A = [ a0; a1,a2,a3,... aK]

= a0 + 1/(a1 + 1/(a2 + 1/(a3 + 1/(... + 1/aK) ...)))


对于0 < p/q < 1,我们将遵循以下算法。

  1. 初始化 Y = 0 = [0] ,Z = 1 = [1] ,k = 0。

  2. Outer loop: The 先决条件 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 = [y0; y1, y2, y3, ... yK] and Z = [y0; y1, y2, y3, ... yK + 1]

    • (- 1) K(Y-X) < 0 < (- 1) K(Z-X) ,或更简单地说,对于 k 偶数,Y < X < Z 和对于 k 奇数,Z < X < Y。

  3. 在不改变数字值的情况下,将连分数的度数延长一步。一般来说,如果最后一项是 yK和 yK + 1,我们将其改为[ ... yK,yK + 1 = & # 8734; ]和[ ... yK,zK + 1 = 1]。现在把 k 增加1。

  4. 内部循环 : 这本质上与@templatetypedef 关于整数的面试问题相同。我们进行两阶段二进制搜索以更接近目标:

  5. 内环1 : yK = & # 8734; ,zK = a,X 在 Y 和 Z 之间。

  6. 双 Z 的最后一项: 计算 M = Z,但使用 mK = 2 * a = 2 * zK

  7. 查询未知数: q = Q (X,M)。

  8. 如果 q = 0,我们就得到了答案,进入第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. 否则,设置 Y = M,然后进入下一步:

  11. 内环2. yk = b,zk = a,X 在 Y 和 Z 之间。

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

  13. 执行二进制搜索: 计算 M,其中 mK = floor ((a + b)/2,查询 q = Q (X,M)。

  14. 如果 q = 0,我们就完成了,进入第17步。

  15. 如果 q 和 Q (X,Y)有相反的符号,这意味着 X 在 Y 和 M 之间,所以设置 Z = M,进入第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. 完成: 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)步骤。

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.

下次有机会的时候,我会发布一个实现这个算法的 python 程序。

编辑 : 还要注意的是,你不必每一步都计算连分数(也就是 O (k) ,有一些部分近似值可以用来计算 O (1)中前一步的下一步)

编辑2 : 部分近似递归定义:

如果 AK = [ a0; a1,a2,a3,... aK] = pK/qK,那么 pK = aKp00 + p01,和 qK = aKq00 + q01。(来源: Niven & Zuckerman,第4版,定理7.3-7.5。另见 06)

Example: [0] = 0/1 = p0/q0, [0; 7] = 1/7 = p1/q1; so [0; 7, 16] = (16*1+0)/(16*7+1) = 16/113 = p2/q2.

这意味着如果除了最后一个项外,两个连分数 Y 和 Z 的项相同,而不包括最后一个项的连分数是 pK-1/qK-1,那么我们可以写 Y = (yKpK-1 + pK-2)/(yKqK-1 + qK-2)和 Z = (zKpK-1 + pK-2)/(zKqK-1 + qK-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 q)而不是 O (log ^ 2q)的大纲:

首先要注意的是,在找到有理数之前,每个新的连分数项的步骤 nk的 # 是 没错2b (a _ k)-1,其中 b (a _ k)是表示 a _ k = ceil (log2(a _ k))所需的位的 # : b (a _ k)步骤扩大二进制搜索的“网”,b (a _ k)-1步骤缩小它)。看看上面的例子,您会注意到步骤的 # 总是1、3、7、15等等。

Now we can use the recurrence relation qK = aKqK-1 + qk-2 and induction to prove the desired result.

让我们这样说: 在达到 kth 项所需的 Nk = sum (nk)步骤之后,q 的值有一个最小值: 对于某些固定常数 A,c,q > = A * 2CN(因此倒过来,我们得到步骤 N 的 # 是 < = (1/c) * log2(q/A) = O (log q))

Base cases:

  • K = 0: q = 1,N = 0,所以 q > = 2N
  • K = 1: 对于 N = 2b-1步,q = a1 > = 2b-1 = 2(N-1)/2 = 2N/2/sqrt (2)。

这意味着 A = 1,c = 1/2可以提供所需的边界。实际上,q 可能是 没有的两倍(反例: [0; 1,1,1,1,1]有一个 phi = (1 + sqrt (5))/2的增长因子) ,所以让我们使用 c = 1/4。

归纳:

  • 对于项 k,qK = aKqk-1 + qK-2。同样,对于本学期所需的 nK = 2b-1步骤,aK > = 2B-1 = 2 < sup > (nK-1)/2 。

    所以 aKqK-1 > = 2 < sup > (NK-1)/2 * qK-1 > = 2 < sup > (nK-1)/2 * A * 2 < sup > NK-1/4 = A * 2 < sup > NK/4 /sqrt (2) * 2 < sup > nK/4 。

这里比较困难的部分是,如果 aK = 1,q 在这一项中可能不会增加很多,我们需要使用 qK-2,但它可能比 qK-1小得多。