你得到一个整数n,你需要在斯特恩的双原子序列中找到它的第一次出现的索引.
序列定义如下:
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
我们可以在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其右子节点的节点来扩充Calkin-Wilf树是很好的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.
现在我们需要找到一个分子小于其分母的节点.那当然是其父母的左子女.那么父母b/q就是b/(q-b).如果
q = c*b + d, 1 <= d < b
我们必须c从节点b/d到左边的孩子时间到达b/q.
等等.
我们可以找到从root(1/1)到p/q节点的方式,使用连续分数(我只考虑这里的简单连续分数)扩展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_11 开头的二进制表示,后跟c_2零,然后c_3是1等,以...结尾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是奇数,k是a[2*j+1]/a[2*j+2] = (a[j] + a[j+1])/a[j+1],因此它是一个正确的孩子.所以,最小的指数k与a[k] = n对应于所有与分子节点的最短路径的最左边的结局n.
最短路径对应于连续分数的扩展n/m,其中0 < m <= n与n[部分必须减少]相互作用,其中部分商的最小总和.
我们需要什么样的长度？由于连分数p/q = [a_0, a_1, ..., a_r]与a_0 > 0和sum
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,我们假设).
通过以下方式调整找到的连续分数扩展[数量不大]:
[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)