Первое появление в диатомической секвенции Стерна

Question

Вы получаете целое число 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

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)

Ключ к этому - дерево Калкин-Уилф.

Начиная с дроби 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) - это последовательность числителей дробей в дереве Калкина-Уилфа, когда они проходят уровень за уровнем, каждый уровень слева направо.

Если мы посмотрим на дерево индексов

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

мы можем легко проверить, что узел по индексу k в дереве Калкин-Уилф несет фракцию 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] снова по определяющим уравнениям.

Все положительные редуцированные дроби встречаются ровно один раз в дереве Калкина-Уилфа (оставлено в качестве упражнения для читателя), поэтому все положительные целые числа встречаются в двухатомной последовательности.

Мы можем найти узел в дереве Калкина-Уилфа по индексу, следуя двоичному представлению индекса, от старшего значащего к младшему, для 1-битного мы переходим к правому дочернему элементу и для 0-битного к левый. (Для этого неплохо дополнить дерево Калкина-Уилфа узлом 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 раз.

Теперь нам нужно найти узел, числитель которого меньше, чем его знаменатель. Это, конечно, левый потомок своего родителя. Родитель 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 ровно, то иди к нужному ребенку 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_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], Длина пути

1. количество битов, необходимых для представления k, а также
2. сумма частных отношений в продолженном разложении дроби.

Теперь, чтобы найти индекс первого появления n > 0 в двухатомной последовательности мы сначала заметим, что наименьший индекс обязательно должен быть нечетным, поскольку a[k] = a[k/2] даже для k, Пусть наименьший индекс будет k = 2*j+1, затем

1. длина RLE k странно,
2. дробь в узле с индексом 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 и сумма

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

1. за 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 а также 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)

Answer 2

Я бы порекомендовал вам прочитать это письмо от Дейкстры, которое объясняет альтернативный способ вычисления этой функции с помощью:

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.

Следовательно, индекс первого появления определенного значения для b может быть вычислен путем нахождения наименьшего n, для которого эта итерация приведет к этому значению b.

Один из способов найти это наименьшее n - это использовать поиск A*, где стоимость - это значение n. Эффективность алгоритма будет зависеть от выбранного вами эвристического алгоритма.

Для эвристики я бы рекомендовал отметить, что:

1. конечное значение всегда будет кратно gcd(a,b) (это может быть использовано для исключения некоторых узлов, которые никогда не могут создать цель)
2. б всегда увеличивается
3. существует максимальная (экспоненциальная) скорость, с которой 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)]