Быстрое вычисление Фибоначчи

python algorithm performance fibonacci

1933 просмотра

4 ответа

2012 Репутация автора

Несколько недель назад я увидел комментарий в Google+, в котором кто-то продемонстрировал прямое вычисление чисел Фибоначчи, которое не было основано на рекурсии и не использовало запоминание. Он просто запомнил последние 2 цифры и продолжал добавлять их. Это алгоритм O (n), но он реализовал его очень чисто. Поэтому я быстро указал, что более быстрый способ - воспользоваться тем, что они могут быть вычислены как степени матрицы [[0,1], [1,1]], и для этого требуется только O (log (N)) вычисление.

Проблема, конечно, в том, что это далеко не оптимально за определенный момент. Это эффективно до тех пор, пока числа не слишком велики, но их длина увеличивается со скоростью N * log (phi) / log (10), где N - это N-е число Фибоначчи, а фи - золотое сечение ((1 + sqrt (5)) / 2 ~ 1,6). Как оказалось, log (phi) / log (10) очень близок к 1/5. Таким образом, можно ожидать, что число N-й Фибоначчи будет иметь примерно N / 5 цифр.

Матричное умножение, даже умножение четных чисел, становится очень медленным, когда числа начинают иметь миллионы или миллиарды цифр. Таким образом, для вычисления F (100 000) потребовалось около 0,03 секунды (в Python), а для F (1000 000) - примерно 5 секунд. Это вряд ли O (log (N)) роста. Моя оценка состояла в том, что этот метод, без улучшений, оптимизирует только вычисления для O ((log (N)) ^ (2.5)) или около того.

Вычисление миллиардного числа Фибоначчи с такой скоростью было бы чрезмерно медленным (даже если бы оно имело всего ~ 1 000 000 000/5 цифр, поэтому оно легко помещалось в 32-разрядную память).

Кто-нибудь знает о реализации или алгоритме, который позволил бы более быстрые вычисления? Возможно, что-то, что позволило бы вычислить триллионное число Фибоначчи.

И просто чтобы быть ясным, я не ищу приближения. Я ищу точное вычисление (до последней цифры).

Редактировать 1: я добавляю код Python, чтобы показать то, что я считаю алгоритмом O ((log N) ^ 2.5)).

from operator import mul as mul
from time import clock

class TwoByTwoMatrix:
    __slots__ = "rows"

    def __init__(self, m):
        self.rows = m

    def __imul__(self, other):
        self.rows = [[sum(map(mul, my_row, oth_col)) for oth_col in zip(*other.rows)] for my_row in self.rows]
        return self

    def intpow(self, i):
        i = int(i)
        result = TwoByTwoMatrix([[long(1),long(0)],[long(0),long(1)]])
        if i <= 0:
            return result
        k = 0
        while i % 2 == 0:
            k +=1
            i >>= 1
        multiplier = TwoByTwoMatrix(self.rows)
        while i > 0:
            if i & 1:
                result *= multiplier
            multiplier *= multiplier # square it
            i >>= 1
        for j in xrange(k):
            result *= result
        return result


m = TwoByTwoMatrix([[0,1],[1,1]])

t1 = clock()
print len(str(m.intpow(100000).rows[1][1]))
t2 = clock()
print t2 - t1

t1 = clock()
print len(str(m.intpow(1000000).rows[1][1]))
t2 = clock()
print t2 - t1

Изменить 2: Похоже, я не учел тот факт, len(str(...))что внесет существенный вклад в общее время выполнения теста. Изменение тестов на

from math import log as log

t1 = clock()
print log(m.intpow(100000).rows[1][1])/log(10)
t2 = clock()
print t2 - t1

t1 = clock()
print log(m.intpow(1000000).rows[1][1])/log(10)
t2 = clock()
print t2 - t1

сократили время выполнения до 0,008 секунды и 0,31 секунды (с 0,03 секунды и 5 секунд при len(str(...))использовании).

Поскольку M = [[0,1], [1,1]] возводится в степень N, это [[F (N-2), F (N-1)], [F (N-1), F (N) ]], другим очевидным источником неэффективности было вычисление (0,1) и (1,0) элементов матрицы, как если бы они были различны. Это (и я переключился на Python3, но Python2.7 раза похожи):

class SymTwoByTwoMatrix():
    # elments (0,0), (0,1), (1,1) of a symmetric 2x2 matrix are a, b, c.
    # b is also the (1,0) element because the matrix is symmetric

    def __init__(self, a, b, c):
        self.a = a
        self.b = b
        self.c = c

    def __imul__(self, other):
        # this multiplication does work correctly because we 
        # are multiplying powers of the same symmetric matrix
        self.a, self.b, self.c = \
            self.a * other.a + self.b * other.b, \
            self.a * other.b + self.b * other.c, \
            self.b * other.b + self.c * other.c
        return self

    def intpow(self, i):
        i = int(i)
        result = SymTwoByTwoMatrix(1, 0, 1)
        if i <= 0:
            return result
        k = 0
        while i % 2 == 0:
            k +=1
            i >>= 1
        multiplier = SymTwoByTwoMatrix(self.a, self.b, self.c)
        while i > 0:
            if i & 1:
                result *= multiplier
            multiplier *= multiplier # square it
            i >>= 1
        for j in range(k):
            result *= result
        return result

рассчитано F (100 000) в 0,006, F (1 000 000) в 0,235 и F (10 000 000) за 9,51 секунды.

Что и следовало ожидать. Он дает результаты на 45% быстрее для самого быстрого теста, и ожидается, что усиление должно асимптотически приближаться к фи / (1 + 2 * фи + фи * фи) ~ 23,6%.

Элемент (0,0) в M ^ N на самом деле является N-вторым числом Фибоначчи:

for i in range(15):
    x = m.intpow(i)
    print([x.a,x.b,x.c])

дает

[1, 0, 1]
[0, 1, 1]
[1, 1, 2]
[1, 2, 3]
[2, 3, 5]
[3, 5, 8]
[5, 8, 13]
[8, 13, 21]
[13, 21, 34]
[21, 34, 55]
[34, 55, 89]
[55, 89, 144]
[89, 144, 233]
[144, 233, 377]
[233, 377, 610]

Я ожидаю, что отсутствие необходимости вычислять элемент (0,0) приведет к увеличению скорости на 1 / (1 + фи + фи * фи) ~ 19%. Но решение lru_cacheF (2N) и F (2N-1), данное Эли Корвиго ниже, на самом деле дает увеличение скорости в 4 раза (то есть на 75%). Таким образом, хотя я не разработал формальное объяснение, я склонен думать, что он кэширует промежутки 1 в двоичном расширении N и выполняет минимальное количество необходимых умножений. Это избавляет от необходимости находить эти диапазоны, предварительно вычислять их, а затем умножать их в нужной точке при расширении N. lru_cacheпозволяет вычислять сверху вниз то, что было бы более сложным вычислением снизу вверх.

И SymTwoByTwoMatrixдля lru_cache-of-F (2N) -and-F (2N-1) требуется примерно в 40 раз больше времени для вычисления каждый раз, когда N увеличивается в 10 раз. Я думаю, что это возможно из-за реализации умножения длинных целых чисел в Python. Я думаю, что умножение больших чисел и их сложение должны быть распараллеливаемыми. Таким образом, многопоточное решение sub-O (N) должно быть возможным, даже если (как утверждает Даниэль Фишер в комментариях) решение F (N) есть Theta(n).

Автор: Dmitry Rubanovich Источник Размещён: 18.07.2016 08:01

Ответы (4)


0 плюса

4909 Репутация автора

Из Википедии ,

Для всех n ≥ 0 число Fn является ближайшим целым числом к ​​phi ^ n / sqrt (5) где phi - золотое сечение. Поэтому его можно найти округлением, то есть использованием функции ближайшего целого

Автор: Aaron Размещён: 18.07.2016 08:13

5 плюса

6682 Репутация автора

Поскольку последовательность Фибоначчи является линейным рецидивом, ее члены могут оцениваться в замкнутой форме. Это включает в себя вычисление мощности, которое может быть выполнено в O (logn) аналогично решению умножения матриц, но постоянные издержки должны быть ниже. Это самый быстрый алгоритм, который я знаю.

выдумка

РЕДАКТИРОВАТЬ

Извините, я пропустил "точную" часть. Другая точная O (log (n)) альтернатива для умножения матриц может быть вычислена следующим образом

fib2

from functools import lru_cache

@lru_cache(None)
def fib(n):
    if n in (0, 1):
        return 1
    if n & 1:  # if n is odd, it's faster than checking with modulo
        return fib((n+1)//2 - 1) * (2*fib((n+1)//2) - fib((n+1)//2 - 1))
    a, b = fib(n//2 - 1), fib(n//2)
    return a**2 + b**2

Это основано на выводе из заметки профессора Эдсгера Дейкстры. Решение использует тот факт, что для вычисления F (2N) и F (2N-1) вам нужно знать только F (N) и F (N-1). Тем не менее, вы все еще имеете дело с арифметикой длинных чисел, хотя издержки должны быть меньше, чем у матричного решения. В Python лучше переписать это в императивном стиле из-за медленной запоминания и рекурсии, хотя я написал это таким образом для ясности функциональной формулировки.

Автор: Eli Korvigo Размещён: 18.07.2016 08:19

0 плюса

8618 Репутация автора

Это слишком долго для комментария, поэтому я оставлю ответ.

Ответ Аарона правильный, и я проголосовал за него, как и вы. Я предоставлю тот же ответ и объясню, почему он не только правильный, но и лучший ответ, опубликованный до сих пор. Формула, которую мы обсуждаем:

формула

Вычисление Φ - это O(M(n)), где M(n)сложность умножения (в настоящее время немного более линейная) и nчисло бит.

Тогда есть степенная функция, которая может быть выражена как log ( O(M(n)•log(n)), multiply ( O(M(n))) и exp ( O(M(n)•log(n)).

Тогда есть квадратный корень ( O(M(n))), деление ( O(M(n))) и последний раунд ( O(n)).

Это делает этот ответ что - то вроде O(n•log^2(n)•log(log(n)))для nбит.


Я не полностью проанализировал алгоритм деления, но если я правильно понял, каждый бит может нуждаться в рекурсии (вам нужно разделить число log(2^n)=nраз), и каждая рекурсия нуждается в умножении. Поэтому это не может быть лучше O(M(n)•n), и это экспоненциально хуже .

Автор: imallett Размещён: 21.07.2016 07:56

1 плюс

2554 Репутация автора

Используя странное уравнение квадратного корня в другом ответе, закрытом от Фибоначчи, вы МОЖЕТЕ точно вычислить k-е число Фибоначчи. Это потому, что $ \ sqrt (5) $ выпадает в конце. Вы просто должны организовать свое умножение, чтобы отслеживать его в то же время.

def rootiply(a1,b1,a2,b2,c):
    ''' multipy a1+b1*sqrt(c) and a2+b2*sqrt(c)... return a,b'''
    return a1*a2 + b1*b2*c, a1*b2 + a2*b1

def rootipower(a,b,c,n):
    ''' raise a + b * sqrt(c) to the nth power... returns the new a,b and c of the result in the same format'''
    ar,br = 1,0
    while n != 0:
        if n%2:
            ar,br = rootiply(ar,br,a,b,c)
        a,b = rootiply(a,b,a,b,c)
        n /= 2
    return ar,br

def fib(k):
    ''' the kth fibonacci number'''
    a1,b1 = rootipower(1,1,5,k)
    a2,b2 = rootipower(1,-1,5,k)
    a = a1-a2
    b = b1-b2
    a,b = rootiply(0,1,a,b,5)
    # b should be 0!
    assert b == 0
    return a/2**k/5

if __name__ == "__main__":
    assert rootipower(1,2,3,3) == (37,30) # 1+2sqrt(3) **3 => 13 + 4sqrt(3) => 39 + 30sqrt(3)
    assert fib(10)==55
Автор: Scott Размещён: 20.01.2017 07:34
Вопросы из категории :
32x32