algorithm - First appearance in Stern's Diatomic Sequance -

you integer n , need find index of first appearance in stern's diatomic sequence.

the sequence defined this:

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

see mathworld.

because n can 400000, it's not idea brute-force it, since time limit 4000 ms.

the sequence pretty odd: first occurrence of 8 21, first occurrence of 6 33.

any ideas how solve this?

maybe might help: oeis

we can solve first occurrence of number in range of 400000 in under 4 seconds:

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) 

the key calkin-wilf tree.

starting fraction 1/1, built rule node fraction a/b, left child carries fraction a/(a+b), , right child fraction (a+b)/b.

                        1/1                        /   \                       /     \                      /       \                   1/2         2/1                  /   \       /   \                1/3   3/2   2/3   3/1 

etc. diatomic sequence (starting @ index 1) sequence of numerators of fractions in calkin-wilf tree, when traversed level level, each level left right.

if @ tree of indices

                         1                         / \                        /   \                       /     \                      2       3                     / \     / \                    4   5   6   7                   / \                   8   9 ... 

we can verify node @ index k in calkin-wilf tree carries fraction a[k]/a[k+1] induction.

that true k = 1 (a[1] = a[2] = 1), , on,

• for k = 2*j have left child of node index j, fraction a[j]/(a[j]+a[j+1]) , a[k] = a[j] , a[k+1] = a[j] + a[j+1] defining equations of sequence.

• for k = 2*j+1 have right child of node index j, fraction (a[j]+a[j+1])/a[j+1] , a[k]/a[k+1] again defining equations.

all positive reduced fractions occur once in calkin-wilf tree (left exercise reader), hence positive integers occur in diatomic sequence.

we can find node in calkin-wilf tree index following binary representation of index, significant bit least, 1-bit go right child , 0-bit left. (for that, nice augment calkin-wilf tree node 0/1 right child 1/1 node, need have step significant set bit of index.)

now, doesn't yet solve problem @ hand.

but, let first solve related problem: reduced fraction p/q, determine index.

suppose p > q. know p/q right child, , parent (p-q)/q. if p-q > q, have again right child, parent (p - 2*q)/q. continuing, if

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

then reach p/q node b/q node going right child a times.

now need find node numerator smaller denominator. of course left child of parent. parent of b/q b/(q-b) then. if

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

we have go left child c times node b/d reach b/q.

and on.

we can find way root (1/1) p/q node using continued fraction (i consider simple continued fractions here) expansion of p/q. let p > q and

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

the continued fraction expansion of p/q ending in 1.

• if r even, go right child a_r times, left a_(r-1) times, right child ... a_1 times left child, , a_0 times right.
• if r odd, first go left child a_r times, a_(r-1) times right ... a_1 times left child, , a_0 times right.

for p < q, must end going left, hence start going left r , start going right odd r.

we have found close connection between binary representation of index , continued fraction expansion of fraction carried node via path root node.

let run-length-encoding of index k be

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

i.e. binary representation of k starts c_1 ones, followed c_2 zeros, c_3 ones etc., , ending c_j

• ones, if k odd - hence j odd;
• zeros, if k - hence j even.

then [c_j, c_(j-1), ..., c_2, c_1] continued fraction expansion of a[k]/a[k+1] length has same parity k (every rational has 2 continued fraction expansions, 1 odd length, other length).

the rle gives path 0/1 node above 1/1 a[k]/a[k+1]. length of path is

1. the number of bits necessary represent k, and
2. the sum of partial quotients in continued fraction expansion.

now, find index of first occurrence of n > 0 in diatomic sequence, first observe smallest index must odd, since a[k] = a[k/2] k. let smallest index k = 2*j+1. then

1. the length of rle of k odd,
2. the fraction @ node index k a[2*j+1]/a[2*j+2] = (a[j] + a[j+1])/a[j+1], hence right child.

so smallest index k a[k] = n corresponds left-most ending of shortest paths node numerator n.

the shortest paths correspond continued fraction expansions of n/m, 0 < m <= n coprime n [the fraction must reduced] smallest sum of partial quotients.

what kind of length need expect? given continued fraction p/q = [a_0, a_1, ..., a_r] a_0 > 0 , sum

s = a_0 + ... + a_r 

the numerator p bounded f(s+1) , denominator q f(s), f(j) j-th fibonacci number. bounds sharp, a_0 = a_1 = ... = a_r = 1 fraction f(s+1)/f(s).

so if f(t) < n <= f(t+1), sum of partial quotients of continued fraction expansion (either of two) >= t. there m such sum of partial quotients of continued fraction expansion of n/m t, not always:

f(5) = 5 < 6 <= f(6) = 8 

and continued fraction expansions of 2 reduced fractions 6/m 0 < m <= 6 are

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

with sum of partial quotients 6. however, smallest possible sum of partial quotients never larger (the largest i'm aware of t+2).

the continued fraction expansions of n/m , n/(n-m) closely related. let's assume m < n/2, , let

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

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

the continued fraction expansion of n/(n-m) is

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

in particular, sum of partial quotients same both.

unfortunately, i'm not aware of way find m smallest sum of partial quotients without brute force, algorithm (i assume n > 2

1. for 0 < m < n/2 coprime n, find continued fraction expansion of n/m, collecting ones smallest sum of partial quotients (the usual algorithm produces expansions last partial quotient > 1, assume that).

2. adjust found continued fraction expansions [those not large in number] following way:

   • if cf [a_0, a_1, ..., a_r] has length, convert [a_0, a_1, ..., a_(r-1), a_r - 1, 1]
   • otherwise, use [1, a_0 - 1, a_1, ..., a_(r-1), a_r - 1, 1]

   (that chooses 1 between n/m , n/(n-m) leading smaller index)

3. reverse continued fractions obtain run-length-encodings of corresponding indices

4. choose smallest among them.

in step 1, useful use smallest sum found far short-cut.

code (haskell, since that's easiest):

module diatomic (diatomic, firstdiatomic, fuscs)  import data.list  strip :: int -> int -> int strip p = go       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               go 1 _ = []         go m p             | m < p*p   = [m]             | r == 0    = p : go (strip p q) (p+2)             | otherwise = go m (p+2)                               (q,r) = m `quotrem` p  contfraclim :: int -> int -> int -> maybe [int] contfraclim = go []       go acc lim n d = case n `quotrem` d of                        (a,b) | lim < -> nothing                              | b == 0  -> (a:acc)                              | otherwise -> go (a:acc) (lim - a) d b  fixupcf :: [int] -> [int] fixupcf [a]     | < 3     = [a]     | otherwise = [1,a-2,1] fixupcf xs     | (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) = : fixend bs fixend _ = error "shouldn't have called fixend empty list"  cfcompare :: [int] -> [int] -> ordering cfcompare (a:bs) (c:ds) = case compare 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                           (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       pfs = primefactors n     (step,ops) = case pfs of                    (2:xs) -> (2,xs)                    _      -> (1,pfs)     start0 = (n-1) `quot` 2     start | n && start0 = start0 - 1           | otherwise             = start0     eligible k = ((/= 0) . (k `rem`)) ops     estimate = length (takewhile (<= fromintegral n) fibs) + 2     check candidates lim k         | k < 1 || n `quot` k >= lim = if null candidates                                          check [] (2*lim) start                                          else minimumby cfcompare candidates         | eligible k = case contfraclim lim n k of                          nothing -> check candidates lim (k-step)                          cf -> let s = sum cf                                     in if s < lim                                          check [fixupcf cf] s (k - step)                                          else check (fixupcf cf : candidates) lim (k-step)         | otherwise = check candidates lim (k-step)