首次出现在斯特恩的双原子序列中

你得到一个整数n，你需要找到它在斯特恩双原子序列中第一次出现的索引。

序列定义如下：

a[0]     = 0
a[1]     = 1
a[2*i]   = a[i]
a[2*i+1] = a[i] + a[i+1]

See 数学世界 http://mathworld.wolfram.com/SternsDiatomicSeries.html.

由于 n 最大可达 400000，因此暴力破解它并不是一个好主意，特别是因为时间限制为 4000 毫秒。

这个序列非常奇怪：第一次出现 8 是 21，但第一次出现 6 是 33。

有什么想法如何解决这个问题吗？

也许这可能有帮助：OEIS http://oeis.org/A020946

我们可以在四秒内轻松解决 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），从那时起，

• for k = 2*j我们有带有索引的节点的左子节点j，所以分数是a[j]/(a[j]+a[j+1]) and a[k] = a[j] and a[k+1] = a[j] + a[j+1]是序列的定义方程。

• for 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

然后我们到达p/q节点从b/q通过转到右子节点a times.

现在我们需要找到一个分子小于分母的节点。那当然是其父母的左孩子。的父母b/q is b/(q-b)然后。如果

q = c*b + d, 1 <= d < b

我们必须去找左边的孩子c距离节点的时间b/d达到b/q.

等等。

我们可以从根源上找到路（1/1）到p/q节点使用连分数（我这里只考虑简单的连分数）展开p/q. Let p > q and

p/q = [a_0, a_1, ..., a_r,1]

的连续分数展开式p/q结束于1.

• If r是偶数，则去右边的孩子a_r次，然后向左a_(r-1)次，然后给正确的孩子......然后a_1次到左边的孩子，最后a_0次向右。
• If r是奇数，那么先去左边的孩子a_r次，那么a_(r-1)向右...然后a_1次到左边的孩子，最后a_0次向右。

For p < q，我们必须结束向左移动，因此开始向左移动r并开始向右走奇数r.

因此，我们发现索引的二进制表示与节点通过从根到节点的路径所携带的分数的连续分数展开之间存在紧密的联系。

让索引的游程编码k be

[c_1, c_2, ..., c_j]           (all c_i > 0)

即二进制表示k以。。开始c_1的，然后是c_2零，那么c_3等等，并以c_j

• 那些，如果k是奇数 - 因此j也是奇数；
• 零，如果k是偶数 - 因此j也是均匀的。

Then [c_j, c_(j-1), ..., c_2, c_1]是连续分式展开式a[k]/a[k+1]其长度具有相同的奇偶性k（每个有理数恰好有两个连分数展开式，一个具有奇数长度，另一个具有偶数长度）。

RLE 给出了从0/1上面的节点1/1 to a[k]/a[k+1]。路径的长度是

1. 表示所需的位数k, and
2. 连分数展开式中的部分商之和。

现在，找到第一次出现的索引n > 0在双原子序列中，我们首先观察到最小的索引必然是奇数，因为a[k] = a[k/2]对于甚至k。设最小索引为k = 2*j+1. Then

1. 角色的长度k is odd,
2. 具有索引的节点处的分数k is a[2*j+1]/a[2*j+2] = (a[j] + a[j+1])/a[j+1]，因此它是一个右孩子。

所以最小的索引k with a[k] = n对应于分子节点的所有最短路径的最左端n.

最短路径对应于连续分数展开式n/m, where 0 < m <= n互质于n[必须减少分数]与最小的部分商之和。

我们需要期望什么样的长度？给定一个连分数p/q = [a_0, a_1, ..., a_r] with a_0 > 0 and sum

s = a_0 + ... + a_r

分子p的边界是F(s+1)和分母q by F(s), where F(j) is the j-th 斐波那契数。边界是尖锐的，因为a_0 = a_1 = ... = a_r = 1分数是F(s+1)/F(s).

So if F(t) < n <= F(t+1)，连分数展开式（两者之一）的部分商之和为>= t。经常有一个m使得连分数展开式的部分商之和n/m正是t， 但不总是：

F(5) = 5 < 6 <= F(6) = 8

以及两个约简分数的连分数展开式6/m with 0 < m <= 6 are

6/1 = [6]          (alternatively [5,1])
6/5 = [1,4,1]      (alternatively [1,5])

部分商之和为 6。然而，最小可能的部分商之和永远不会大很多（我知道的最大的是t+2).

的连续分数展开n/m and n/(n-m)密切相关。我们假设m < n/2， 然后让

n/m = [a_0, a_1, ..., a_r]

Then 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) is

n/(n-m) = [1, a_0 - 1, a_1, ..., a_r]

特别是，两者的部分商之和是相同的。

不幸的是，我不知道如何找到m具有最小的部分商之和，无需暴力破解，所以该算法是（我假设n > 2

1. for 0 < m < n/2互质于n，求出连分式展开式n/m，收集部分商之和最小的那些（通常的算法产生的展开式的最后一个部分商是> 1，我们假设）。

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 and 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)