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

Question

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

序列定义如下:

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

见 MathWorld 。

因为 n 可以达到 400000，所以暴力破解它不是一个好主意，特别是因为时间限制是 4000 毫秒。

序列很奇怪:第一次出现 8 是 21，但第一次出现 6 是 33。

任何想法如何解决这个问题？

也许这可能会有所帮助: OEIS

Answer 1

我们可以在四秒内轻松解决第一次出现 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 的节点，其左 child 携带分数 a/(a+b) ，其右 child 携带分数 (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 ...

我们可以很容易地通过归纳验证 Calkin-Wilf 树中索引 k 处的节点携带分数 a[k]/a[k+1]。
k = 1 ( a[1] = a[2] = 1 ) 显然是这样，从那时起，

对于 k = 2*j 我们有索引为 j 的节点的左 child ，所以分数是 a[j]/(a[j]+a[j+1]) 和 a[k] = a[j] 和 a[k+1] = a[j] + a[j+1] 是序列的定义方程。

对于 k = 2*j+1 我们有索引为 j 的节点的右 child ，所以分数是 (a[j]+a[j+1])/a[j+1] 并且通过定义方程再次是 a[k]/a[k+1] 。

所有正归约分数在 Calkin-Wilf 树中只出现一次(留给读者作为练习)，因此所有正整数都出现在双原子序列中。

我们可以按照索引的二进制表示从索引中找到 Calkin-Wilf 树中的节点，从最高有效位到最低有效位，对于 1 位我们去右 child ，对于 0 位到左边。 (为此，最好用节点 0/1 扩充 Calkin-Wilf 树，该节点的右子节点是 1/1 节点，这样我们就需要对索引的最重要的设置位有一个步骤。)

现在，这对解决手头的问题还没有多大帮助。

但是，让我们首先解决一个相关的问题:对于约简分数 p/q ，确定其索引。

假设 p > q 。那么我们知道 p/q 是一个右 child ，它的 parent 是 (p-q)/q 。如果也是 p-q > q ，我们又有一个右 child ，其 parent 是 (p - 2*q)/q 。继续，如果

p = a*q + b, 1 <= b < q

然后我们通过转到右子节点 p/q 次，从 b/q 节点到达 a 节点。

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

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

我们必须从节点 c 到左子节点 b/d 次才能到达 b/q 。

等等。

我们可以使用 1/1 的连分数(这里只考虑简单连分数)展开式找到从根节点 ( p/q ) 到 p/q 节点的路径。让 p > q 和

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

p/q 的连分数展开式以 1 结尾。

如果 r 是偶数，则到右 child a_r 次，然后到左 a_(r-1) 次，然后到右 child ... 然后 a_1 次到左 child ，最后a_0 次向右。

如果 r 是奇数，那么首先去左 child a_r 次，然后 a_(r-1) 次到右边......然后 a_1 次到左 child ，最后 a_0 次到正确的。

对于 p < q ，我们必须结束向左走，因此开始向左走偶数 r 并开始向右走奇数 r 。

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

让索引 k的游程编码为

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

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

个，如果 k 是奇数 - 因此 j 也是奇数；

为零，如果 k 是偶数 - 因此 j 也是偶数。

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

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

表示 k 和

所需的位数

连分数展开式中部分商的总和。

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

k 的 RLE 的长度是奇数，

索引为 k 的节点处的分数是 a[2*j+1]/a[2*j+2] = (a[j] + a[j+1])/a[j+1] ，因此它是右 child 。

因此，具有 k 的最小索引 a[k] = n 对应于所有到具有分子 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 -th Fibonacci 数。界限是尖锐的，对于 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) 之间选择一个导致较小的索引)

反转连分数以获得相应索引的行程编码

选择其中最小的。

在第 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)