2个回答
8
我们可以在不到四秒钟的时间内轻松地解决400000范围内第一次出现的数字：
关键是Calkin-Wilf树。
从分数1/1开始，它是按照以下规则构建的：对于一个带有分数a/b的节点，它的左子节点为分数a/(a+b)，右子节点为分数(a+b)/b。
等等。您的文本似乎是计算机科学中的技术性术语。您能否提供更多上下文，以便我更好地翻译它们？
我们可以通过归纳法轻松验证 Calkin-Wilf 树中索引为 k 的节点携带分数 a[k]/a[k+1]。显然对于 k=1（a[1]=a[2]=1），这是正确的，从那时起，
- 对于 k = 2*j，我们有索引为 j 的节点的左子节点，因此分数为 a[j]/(a[j]+a[j+1])，并且 a[k]=a[j] 和 a[k+1]=a[j]+a[j+1] 是该序列的定义方程。 - 对于 k = 2*j+1，我们有索引为 j 的节点的右子节点，因此分数为 (a[j]+a[j+1])/a[j+1]，根据定义方程，这仍然是 a[k]/a[k+1]。
所有正约分分数在 Calkin-Wilf 树中恰好出现一次（留给读者作为练习），因此所有正整数都出现在二元序列中。
我们可以通过按照索引的二进制表示从最高位到最低位进行遍历来找到 Calkin-Wilf 树中的节点，对于一个 1 位，我们转到右子节点，对于一个 0 位，我们转到左子节点。（为此，将 Calkin-Wilf 树与节点 0/1 扩充，其右子节点是 1/1 节点，因此我们需要在索引的最高位设置一个步骤。）
现在，这还不能很好地解决手头的问题。
但是，让我们先解决一个相关问题：对于约分分数 p/q，确定其索引。
假设 p>q。那么我们知道 p/q 是右子节点，其父节点是 (p-q)/q。如果 p-q > q，我们再次有一个右子节点，其父节点为 (p-2*q)/q。继续，如果
然后我们通过向右子节点移动
现在我们需要找到一个分子小于分母的节点。这当然是其父节点的左子节点。
我们需要从节点b/d向左子节点移动次，才能到达b/q节点。
即，二进制表示为
分数的分子p被F（s+1）限制，分母q被F（s）限制，其中F（j）是第j个斐波那契数。 对于a0 = a1 = ... = ar = 1，分数是F（s +1）/ F（s）。因此，如果F（t）< n <= F（t + 1），则连分数展开的部分商之和（任一方）> = t。通常有一个m，使得n / m的连分数展开的部分商之和恰好为t，但并非总是如此。
两个约分分数的连分数展开式为
“求和为6的偏商之和。然而，偏商之和的最小可能值不会比这个大很多（我知道的最大值是 t+2）。
数值 n/m 和 n/(n-m) 的连分数展开式密切相关。我们假设 m < n/2，并且让”
然后，
"
特别地，两者的部分商之和相同。
不幸的是，我不知道有没有一种方法可以在不使用暴力的情况下找到部分商之和最小的m，所以算法是（我假设n>2）：
1. 对于01的展开式，我们假设这个条件）。 2. 调整找到的连分数展开式[它们的数量不多]如下： - 如果CF [a_0, a_1, ..., a_r] 的长度为偶数，则将其转换为 [a_0, a_1, ..., a_(r-1), a_r - 1, 1] - 否则，使用 [1, a_0 - 1, a_1, ..., a_(r-1), a_r - 1, 1] （选择n/m和n/(n-m)中导致更小索引的那一个） 3. 反转连分数以获得相应索引的运行长度编码 4. 在它们中选择最小值。
在步骤1中，使用迄今为止找到的最小和来进行快捷操作非常有用。
代码（Haskell，因为它最容易）：
Prelude Diatomic> firstDiatomic 400000
363490989
(0.03 secs, 26265328 bytes)
Prelude Diatomic> map firstDiatomic [400000 .. 400100]
[363490989,323659475,580472163,362981813,349334091,355685483,346478235,355707595
,291165867,346344083,347155797,316314293,576398643,315265835,313171245,355183267
,315444051,315970205,575509833,311741035,340569429,313223987,565355925,296441165
,361911645,312104147,557145429,317106853,323637939,324425077,610613547,311579309
,316037811,311744107,342436533,348992869,313382235,325406123,355818699,312128723
,347230875,324752171,313178421,312841811,313215645,321754459,576114987,325793195
,313148763,558545581,355294101,359224397,345462093,307583675,355677549,312120731
,341404245,316298389,581506779,345401947,312109779,316315061,315987123,313447771
,361540179,313878107,304788843,325765547,316036275,313731751,355635795,312035947
,346756533,313873883,349358379,357393763,559244877,313317739,325364139,312128107
,580201947,358182323,314944173,357403987,584291115,312158827,347448723,363246413
,315935571,349386085,315929427,312137323,357247725,313207657,320121429,356954923
,557139285,296392013,576042123,311726765,296408397]
(2.45 secs, 3201358192 bytes)
关键是Calkin-Wilf树。
从分数1/1开始，它是按照以下规则构建的：对于一个带有分数a/b的节点，它的左子节点为分数a/(a+b)，右子节点为分数(a+b)/b。
1/1
/ \
/ \
/ \
1/2 2/1
/ \ / \
1/3 3/2 2/3 3/1
等等。您的文本似乎是计算机科学中的技术性术语。您能否提供更多上下文，以便我更好地翻译它们？
1
/ \
/ \
/ \
2 3
/ \ / \
4 5 6 7
/ \
8 9 ...
我们可以通过归纳法轻松验证 Calkin-Wilf 树中索引为 k 的节点携带分数 a[k]/a[k+1]。显然对于 k=1（a[1]=a[2]=1），这是正确的，从那时起，
- 对于 k = 2*j，我们有索引为 j 的节点的左子节点，因此分数为 a[j]/(a[j]+a[j+1])，并且 a[k]=a[j] 和 a[k+1]=a[j]+a[j+1] 是该序列的定义方程。 - 对于 k = 2*j+1，我们有索引为 j 的节点的右子节点，因此分数为 (a[j]+a[j+1])/a[j+1]，根据定义方程，这仍然是 a[k]/a[k+1]。
所有正约分分数在 Calkin-Wilf 树中恰好出现一次（留给读者作为练习），因此所有正整数都出现在二元序列中。
我们可以通过按照索引的二进制表示从最高位到最低位进行遍历来找到 Calkin-Wilf 树中的节点，对于一个 1 位，我们转到右子节点，对于一个 0 位，我们转到左子节点。（为此，将 Calkin-Wilf 树与节点 0/1 扩充，其右子节点是 1/1 节点，因此我们需要在索引的最高位设置一个步骤。）
现在，这还不能很好地解决手头的问题。
但是，让我们先解决一个相关问题：对于约分分数 p/q，确定其索引。
假设 p>q。那么我们知道 p/q 是右子节点，其父节点是 (p-q)/q。如果 p-q > q，我们再次有一个右子节点，其父节点为 (p-2*q)/q。继续，如果
p = a*q + b, 1 <= b < q
然后我们通过向右子节点移动
a次，从b/q节点到达p/q节点。现在我们需要找到一个分子小于分母的节点。这当然是其父节点的左子节点。
b/q的父节点是b/(q-b)。如果q = c*b + d, 1 <= d < b
我们需要从节点b/d向左子节点移动次，才能到达b/q节点。
以此类推。
我们可以使用p/q的连分数（这里仅考虑简单连分数）展开来找到从根节点（1/1）到p/q节点的路径。假设p>q，且
p/q = [a_0, a_1, ..., a_r,1]
以p/q的连分数展开以1结尾为条件。
- 如果
r是偶数，则先向右子节点走a_r次，然后向左边走a_(r-1)次，再向右子节点...最后向左子节点走a_1次，最终向右子节点走a_0次。 - 如果
r是奇数，则先向左子节点走a_r次，然后向右边走a_(r-1)次...最后向左子节点走a_1次，最终向右子节点走a_0次。
对于p < q，我们必须最后向左边走，因此对于偶数r开始向左边走，对于奇数r开始向右边走。
我们因此找到了索引的二进制表示和通过从根到节点的路径所携带的分数的连分数展开之间的紧密联系。
设索引k的运行长度编码为
[c_1, c_2, ..., c_j] (all c_i > 0)
即，二进制表示为
k以c_1个1开头，后跟c_2个零，然后是c_3个1等，最后以c_j结束。
- 如果
k是奇数，则为1-因此j也是奇数； - 如果
k是偶数，则为0-因此j也是偶数。
然后，[c_j, c_(j-1), ..., c_2, c_1]是a[k]/a[k+1]的连分数展开式，其长度与k具有相同的奇偶性（每个有理数恰好有两个连分数展开式，一个奇数长度，另一个偶数长度）。
RLE给出了从1/1上方的0/1节点到a[k]/a[k+1]的路径。路径的长度为：
- 表示
k所需的位数，以及 - 连分数展开中部分商的总和。
现在，要找到第一次出现n>0的双调序列索引，我们首先观察到最小索引必须是奇数，因为对于偶数k，a[k] = a[k/2]。让最小索引为k=2*j+1。然后
k的RLE长度为奇数，- 索引为
k的节点处的分数为a[2*j+1]/a[2*j+2] = (a[j]+a[j+1])/a[j+1]，因此它是右子节点。
因此，具有a[k]=n的最小索引k对应于所有最短路径的最左端，这些路径指向分子为n的节点。
最短路径对应于n/m的连分数展开式，其中0 < m <= n与n互质[分数必须被约简]，并且其部分商之和最小。
我们需要期望什么样的长度？给定连分数p/q = [a_0, a_1, ..., a_r]，其中a_0 > 0且总和
s = a_0 + ... + a_r
分数的分子p被F（s+1）限制，分母q被F（s）限制，其中F（j）是第j个斐波那契数。 对于a0 = a1 = ... = ar = 1，分数是F（s +1）/ F（s）。因此，如果F（t）< n <= F（t + 1），则连分数展开的部分商之和（任一方）> = t。通常有一个m，使得n / m的连分数展开的部分商之和恰好为t，但并非总是如此。
F(5) = 5 < 6 <= F(6) = 8
两个约分分数的连分数展开式为
6/m，其中0 < m <= 6。6/1 = [6] (alternatively [5,1])
6/5 = [1,4,1] (alternatively [1,5])
“求和为6的偏商之和。然而，偏商之和的最小可能值不会比这个大很多（我知道的最大值是 t+2）。
数值 n/m 和 n/(n-m) 的连分数展开式密切相关。我们假设 m < n/2，并且让”
n/m = [a_0, a_1, ..., a_r]
然后，
a_0 >= 2，(n-m)/m = [a_0 - 1, a_1, ..., a_r]
and since
n/(n-m) = 1 + m/(n-m) = 1 + 1/((n-m)/m)
"
n/(n-m) 的连分数展开式为：
n/(n-m) = [1, a_0 - 1, a_1, ..., a_r]
特别地，两者的部分商之和相同。
不幸的是，我不知道有没有一种方法可以在不使用暴力的情况下找到部分商之和最小的m，所以算法是（我假设n>2）：
1. 对于01的展开式，我们假设这个条件）。 2. 调整找到的连分数展开式[它们的数量不多]如下： - 如果CF [a_0, a_1, ..., a_r] 的长度为偶数，则将其转换为 [a_0, a_1, ..., a_(r-1), a_r - 1, 1] - 否则，使用 [1, a_0 - 1, a_1, ..., a_(r-1), a_r - 1, 1] （选择n/m和n/(n-m)中导致更小索引的那一个） 3. 反转连分数以获得相应索引的运行长度编码 4. 在它们中选择最小值。
在步骤1中，使用迄今为止找到的最小和来进行快捷操作非常有用。
代码（Haskell，因为它最容易）：
module Diatomic (diatomic, firstDiatomic, fuscs) where
import Data.List
strip :: Int -> Int -> Int
strip p = go
where
go n = case n `quotRem` p of
(q,r) | r == 0 -> go q
| otherwise -> n
primeFactors :: Int -> [Int]
primeFactors n
| n < 1 = error "primeFactors: non-positive argument"
| n == 1 = []
| n `rem` 2 == 0 = 2 : go (strip 2 (n `quot` 2)) 3
| otherwise = go n 3
where
go 1 _ = []
go m p
| m < p*p = [m]
| r == 0 = p : go (strip p q) (p+2)
| otherwise = go m (p+2)
where
(q,r) = m `quotRem` p
contFracLim :: Int -> Int -> Int -> Maybe [Int]
contFracLim = go []
where
go acc lim n d = case n `quotRem` d of
(a,b) | lim < a -> Nothing
| b == 0 -> Just (a:acc)
| otherwise -> go (a:acc) (lim - a) d b
fixUpCF :: [Int] -> [Int]
fixUpCF [a]
| a < 3 = [a]
| otherwise = [1,a-2,1]
fixUpCF xs
| even (length xs) = case xs of
(1:_) -> fixEnd xs
(a:bs) -> 1 : (a-1) : bs
| otherwise = case xs of
(1:_) -> xs
(a:bs) -> 1 : fixEnd ((a-1):bs)
fixEnd :: [Int] -> [Int]
fixEnd [a,1] = [a+1]
fixEnd [a] = [a-1,1]
fixEnd (a:bs) = a : fixEnd bs
fixEnd _ = error "Shouldn't have called fixEnd with an empty list"
cfCompare :: [Int] -> [Int] -> Ordering
cfCompare (a:bs) (c:ds) = case compare a c of
EQ -> cfCompare ds bs
cp -> cp
fibs :: [Integer]
fibs = 0 : 1 : zipWith (+) fibs (tail fibs)
toNumber :: [Int] -> Integer
toNumber = foldl' ((+) . (*2)) 0 . concat . (flip (zipWith replicate) $ cycle [1,0])
fuscs :: Integer -> (Integer, Integer)
fuscs 0 = (0,1)
fuscs 1 = (1,1)
fuscs n = case n `quotRem` 2 of
(q,r) -> let (a,b) = fuscs q
in if r == 0
then (a,a+b)
else (a+b,b)
diatomic :: Integer -> Integer
diatomic = fst . fuscs
firstDiatomic :: Int -> Integer
firstDiatomic n
| n < 0 = error "Diatomic sequence has no negative terms"
| n < 2 = fromIntegral n
| n == 2 = 3
| otherwise = toNumber $ bestCF n
bestCF :: Int -> [Int]
bestCF n = check [] estimate start
where
pfs = primeFactors n
(step,ops) = case pfs of
(2:xs) -> (2,xs)
_ -> (1,pfs)
start0 = (n-1) `quot` 2
start | even n && even start0 = start0 - 1
| otherwise = start0
eligible k = all ((/= 0) . (k `rem`)) ops
estimate = length (takeWhile (<= fromIntegral n) fibs) + 2
check candidates lim k
| k < 1 || n `quot` k >= lim = if null candidates
then check [] (2*lim) start
else minimumBy cfCompare candidates
| eligible k = case contFracLim lim n k of
Nothing -> check candidates lim (k-step)
Just cf -> let s = sum cf
in if s < lim
then check [fixUpCF cf] s (k - step)
else check (fixUpCF cf : candidates) lim (k-step)
| otherwise = check candidates lim (k-step)
- Daniel Fischer
6
+1：很好，我认为你的分子和分母与我的b、a相同，但这个答案提供了更清晰的理解和更快的算法。（@dragostis：在我看来，这应该被标记为最佳答案） - Peter de Rivaz
干得好！我很好奇你是否愿意分享一下你得到这些时间的硬件/设置细节。我在1.6Ghz的Windows上运行了这个程序，firstDiatomic 400000 => 2.75秒，map firstDiatomic [400000 .. 400100] => (452.25秒，31745121676字节)；使用fpcomplete.com，400000只需要1.56009秒，而map超时了... - Thell
fpcomplete.com上的地图已完成，用时242.147秒。 - Thell
一个相当标准的2.3 GHz Core i5，运行64位Linux。当然，我已经使用优化编译了模块，解释代码在这里运行也很慢。不过，以1.21 resp 189.4秒的速度并不像你的机器那样慢。 - Daniel Fischer
谢谢 Daniel，现在我对正在使用R和Rcpp的解决方案感到非常满意，它可以获得400000的计时为7ms，[400000..400100]地图的计时为2.4s。 800000..800020花费了2.6秒，仅供娱乐。 - Thell
质数对于处理互质数是致命的。 - Thell
2
我建议你阅读这篇来自Dijkstra的信件，其中解释了通过以下方式计算此函数的另一种方法：
n, a, b := N, 1, 0;
do n ≠ 0 and even(n) → a, n:= a + b, n/2
odd(n) → b, n:= b + a, (n-1)/2
od {b = fusc(N)}
这个算法以a，b = 1,0开头，并且有效地使用了N的连续位（从最低有效位到最高有效位）来增加a和b，最终结果是b的值。
因此，可以通过找到使迭代产生该值b的最小n来计算特定b值的第一次出现的索引。
一种找到这个最小n的方法是使用A*搜索，其中成本是n的值。算法的效率取决于您选择的启发式算法。
对于启发式算法，我建议注意以下几点：
- 最终值始终是gcd（a，b）的倍数（这可以用来排除永远无法产生目标的一些节点）
- b始终在增加
- b可以增加的最大（指数）速率取决于当前a的值
编辑
下面是一些示例Python代码，以说明A *方法。
from heapq import *
def gcd(a,b):
while a:
a,b=b%a,a
return b
def heuristic(node,goal):
"""Estimate least n required to make b==goal"""
n,a,b,k = node
if b==goal: return n
# Otherwise needs to have at least one more bit set
# Improve this heuristic to make the algorithm faster
return n+(1<<k)
def diatomic(goal):
"""Return index of first appearance of n in Stern's Diatomic sequence"""
start=0,1,0,0
f_score=[] # This is used as a heap
heappush(f_score, (0,start) )
while 1:
s,node = heappop(f_score)
n,a,b,k = node
if b==goal:
return n
for node in [ (n,a+b,b,k+1),(n+(1<<k),a,b+a,k+1) ]:
n2,a2,b2,k2 = node
if b2<=goal and (goal%gcd(a2,b2))==0:
heappush(f_score,(heuristic(node,goal),node))
print [diatomic(n) for n in xrange(1,10)]
- Peter de Rivaz
4
你能不能再明确一下，你怎么应用A*算法？从这个表达式中很难看出搜索树。 - dragostis
解决方案似乎是合理的，但我仍然不确定它是否能在400000范围内的4秒内找到任何解决方案。 - dragostis
我对算法符号表示不是很熟悉 --
od 在 od {b = fusc(N)} 中是什么意思？ - גלעד ברקן那就是 do 块的结束。 - dragostis
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接