你得到一个整数 n，你需要找到它在斯特恩双原子序列中第一次出现的索引。
序列定义如下：
a[0] = 0
a[1] = 1
a[2*i] = a[i]
a[2*i+1] = a[i] + a[i+1]
见数学世界。
因为 n 可以达到 400000，所以暴力破解它不是一个好主意，特别是因为时间限制是 4000 毫秒。
这个序列很奇怪：8 的第一次出现是 21，但 6 的第一次出现是 33。
任何想法如何解决这个问题？
也许这可能会有所帮助：OEIS
问问题
977 次
2 回答
8
我们可以在 4 秒内轻松求解 400000 范围内的数字的第一次出现：
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 开始）是 Calkin-Wilf 树中分数的分子序列，当它逐级遍历时，每一级从左到右。
如果我们查看索引树
1
/ \
/ \
/ \
2 3
/ \ / \
4 5 6 7
/ \
8 9 ...
我们可以通过归纳轻松验证kCalkin-Wilf 树中索引处的节点是否携带分数。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 位到左边。（为此，用一个0/1右子节点为该1/1节点的节点来扩充 Calkin-Wilf 树是很好的，这样我们就需要为索引的最高有效位设置一个步骤。）
现在，这对解决手头的问题没有多大帮助。
但是，让我们首先解决一个相关问题：对于约化分数p/q，确定它的索引。
假设p > q. 然后我们知道这p/q是一个正确的孩子，它的父母是(p-q)/q。如果也是p-q > q，我们又有一个正确的孩子，其父母是(p - 2*q)/q。继续，如果
p = a*q + b, 1 <= b < q
然后我们通过去正确的子时间p/q从节点到达节点。b/qa
现在我们需要找到一个分子小于分母的节点。那当然是其父母的左孩子。的父母b/q是b/(q-b)那么。如果
q = c*b + d, 1 <= d < b
我们必须c从节点b/d到左孩子次才能到达b/q。
等等。
我们可以使用 的连分数（我在这里只考虑简单的连分数）展开找到从根 ( 1/1) 到节点的路。让和p/qp/qp > 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)
即二进制表示以1k开头c_1，然后是c_20，然后c_3是 1 等，并以c_j
- 一个，如果
k是奇数 - 因此j也是奇数； - 零，如果
k是偶数 - 因此j也是偶数。
然后[c_j, c_(j-1), ..., c_2, c_1]是a[k]/a[k+1]其长度具有相同奇偶性的连分数展开式k（每个有理数恰好有两个连分数展开式，一个长度为奇数，另一个长度为偶数）。
0/1RLE 给出了从上面的节点1/1到的路径a[k]/a[k+1]。路径的长度是
- 表示 所需的位数
k，以及 - 连分数展开式中的部分商的总和。
现在，要找到n > 0双原子序列中第一次出现 的索引，我们首先观察到最小索引必须是奇数，因为a[k] = a[k/2]对于偶数k。设最小索引为k = 2*j+1。然后
- RLE 的长度
k是奇数， - 具有索引的节点处的分数
k是a[2*j+1]/a[2*j+2] = (a[j] + a[j+1])/a[j+1]，因此它是一个右孩子。
所以最小的索引k对应a[k] = n于所有最短路径的最左端到一个带有 numerator 的节点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-个斐波那契数。界限很明确，因为a_0 = a_1 = ... = a_r = 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]
并且因为
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
对于 与 的
0 < m < n/2互质n，求 的连分数展开式n/m，收集部分商之和最小的那些（通常的算法产生最后一个部分商为 的展开式> 1，我们假设这一点）。通过以下方式调整找到的连分数展开式[数量不大]：
- 如果 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)导致较小的索引）- 如果 CF
反转连分数以获得相应索引的游程编码
选择其中最小的。
在步骤 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)
于 2013-05-05T03:06:09.463 回答
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 的值。
因此，可以通过找到最小的 n 来计算 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)]
于 2013-05-04T19:22:24.730 回答