Cómputo rápido de Fibonacci

Vi un comentario en Google+ hace unas semanas en el que alguien demostró un cálculo directo de los números de Fibonacci que no se basaba en la recursión y no usaba la memoria. Efectivamente solo recordó los últimos 2 números y siguió agregándolos. Este es un algoritmo O (n), pero lo implementó de manera muy limpia. Así que rápidamente señalé que una forma más rápida es aprovechar el hecho de que se pueden calcular como potencias de [[0,1], [1,1]] matriz y solo requiere una O (log (N)) cálculo.

El problema, por supuesto, es que esto está lejos de ser óptimo pasado un cierto punto. Es eficiente siempre que los números no sean demasiado grandes, pero crecen en longitud a la tasa de N * log (phi) / log (10), donde N es el número Nth Fibonacci y phi es la proporción áurea ((1 + sqrt (5)) / 2 ~ 1.6). Como resultado, log (phi) / log (10) está muy cerca de 1/5. Por lo tanto, se puede esperar que el número Nth Fibonacci tenga aproximadamente N / 5 dígitos.

La multiplicación de matrices, heck incluso la multiplicación de números, se vuelve muy lenta cuando los números comienzan a tener millones o billones de dígitos. Así que la F (100,000) tomó cerca de .03 segundos para calcular (en Python), mientras que la F (1000,000) tomó aproximadamente 5 segundos. Esto es apenas un crecimiento O (log (N)). Mi estimación fue que este método, sin mejoras, solo optimiza el cálculo para ser O ((log (N)) ^ (2.5)) o menos.

Cálculo del número de Fibonacci número mil millones, a este ritmo, sería prohibitivamente lento (aunque solo tendría ~ 1,000,000,000 / 5 dígitos, por lo que cabe fácilmente en la memoria de 32 bits).

¿Alguien sabe de una implementación o algoritmo que permita un cálculo más rápido? Tal vez algo que permita el cálculo de un trillón de Fibonacci.

Y para ser claros, no estoy buscando una aproximación. Estoy buscando el cálculo exacto (hasta el último dígito).

Edición 1: estoy agregando el código Python para mostrar lo que creo que es el algoritmo 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 >= 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 

Edición 2: Parece que no tuve en cuenta el hecho de que len(str(...)) haría una contribución significativa al tiempo de ejecución general de la prueba. Cambiando pruebas a

 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 

acortó los tiempos de ejecución a .008 segundos y .31 segundos (de .03 segundos y 5 segundos cuando se usó len(str(...)) ).

Debido a que M = [[0,1], [1,1]] elevado a potencia N es [[F (N-2), F (N-1)], [F (N-1), F (N) ]], la otra fuente obvia de ineficiencia fue el cálculo de (0,1) y (1,0) elementos de la matriz como si fueran distintos. Esto (y cambié a Python3, pero Python2.7 veces son similares):

 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 >= 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 

calculó F (100,000) en .006, F (1,000,000) en .235 y F (10,000,000) en 9.51 segundos.

Lo que es de esperar. Está produciendo resultados 45% más rápido para la prueba más rápida y se espera que la ganancia se aproxime asintóticamente phi / (1 + 2 * phi + phi * phi) ~ 23.6%.

El elemento (0,0) de M ^ N es en realidad el número N-2nd de Fibonacci:

 for i in range(15): x = m.intpow(i) print([xa,xb,xc]) 

da

 [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] 

Espero que no tener que calcular el elemento (0,0) produzca una aceleración adicional de 1 / (1 + phi + phi * phi) ~ 19%. Pero la lru_cache de F (2N) y F (2N-1) dada por Eli Korvigo a continuación en realidad da una velocidad de hasta 4 veces (es decir, 75%). Entonces, si bien no he elaborado una explicación formal, me siento tentado a pensar que almacena los intervalos de 1 dentro de la expansión binaria de N y hace el número mínimo de multiplicaciones necesarias. Lo que evita la necesidad de encontrar esos rangos, precomputarlos y luego multiplicarlos en el punto correcto en la expansión de N. lru_cache permite un cálculo de arriba a abajo de lo que habría sido un cómputo más complicado de abajo a arriba.

Tanto SymTwoByTwoMatrix como lru_cache-of-F (2N) -y-F (2N-1) tardan aproximadamente 40 veces más en calcular cada vez que N crece 10 veces. Creo que posiblemente se deba a la implementación de Python de la multiplicación de las largas intenciones. Creo que la multiplicación de grandes números y su sum deberían ser paralelizables. Por lo tanto, una solución sub-O (N) de múltiples subprocesos debería ser posible aunque (como indica Daniel Fisher en los comentarios) la solución F (N) es Theta(n) .

Dado que la secuencia de Fibonacci es una recurrencia lineal, sus miembros pueden evaluarse en forma cerrada. Esto implica calcular una potencia, que se puede hacer en O (logn) de manera similar a la solución de multiplicación de matrices, pero la sobrecarga constante debe ser menor. Ese es el algoritmo más rápido que conozco.

mentira

EDITAR

Lo siento, me perdí la parte “exacta”. Otra alternativa exacta de O (log (n)) para la multiplicación de matrices se puede calcular de la siguiente manera

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 

Esto se basa en la derivación de una nota del profesor Edsger Dijkstra. La solución explota el hecho de que para calcular F (2N) y F (2N-1) solo necesita conocer F (N) y F (N-1). Sin embargo, usted todavía está tratando con aritméticos de números largos, aunque la sobrecarga debe ser más pequeña que la de la solución basada en matriz. En Python, es mejor que vuelvas a escribir esto en un estilo imperativo debido a la lenta memorización y la recursión, aunque lo escribí de esta manera para la claridad de la formulación funcional.

Usando la ecuación de raíz cuadrada extraña en la otra respuesta, forma cerrada, FIB, usted puede calcular exactamente el número de fibonacci kth. Esto se debe a que $ \ sqrt (5) $ se cae al final. Solo tienes que organizar tu multiplicación para seguirla mientras tanto.

 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 

De Wikipedia ,

Para todos n ≥ 0, el número Fn es el número entero más cercano a phi ^ n / sqrt (5) donde phi es la proporción áurea. Por lo tanto, se puede encontrar redondeando, es decir, mediante el uso de la función de entero más cercana

Esto es demasiado largo para un comentario, así que dejaré una respuesta.

La respuesta de Aaron es correcta, y la he votado, como tú deberías. Proporcionaré la misma respuesta y explicaré por qué no solo es correcta, sino también la mejor respuesta publicada hasta ahora. La fórmula que estamos discutiendo es:

fórmula

El cálculo de Φ es O(M(n)) , donde M(n) es la complejidad de la multiplicación ( actualmente un poco más que la lineal) y n es el número de bits.

Luego hay una función de poder, que se puede express como un registro ( O(M(n)•log(n) ), un multiplicador ( O(M(n)) ), y un exp ( O(M(n)•log(n) ).

Luego hay una raíz cuadrada ( O(M(n)) ), una división ( O(M(n)) ) y una ronda final ( O(n) ).

Esto hace que esta respuesta sea algo como O(n•log^2(n)•log(log(n))) para n bits.


No he analizado a fondo el algoritmo de división, pero si estoy leyendo esto correctamente, cada bit puede necesitar una recursión (es necesario dividir el número de log(2^n)=n veces) y cada recursión necesita una multiplicación. Por lo tanto, no puede ser mejor que O(M(n)•n) , y eso es exponencialmente peor .