performance - serie - enésimo número de fibonacci en tiempo sublineal
sucesion de fibonacci conejos (13)
Aquí está mi versión recursiva que recurses log (n) veces. Creo que es más fácil de leer en la forma recursiva:
def my_fib(x):
if x < 2:
return x
else:
return my_fib_helper(x)[0]
def my_fib_helper(x):
if x == 1:
return (1, 0)
if x % 2 == 1:
(p,q) = my_fib_helper(x-1)
return (p+q,p)
else:
(p,q) = my_fib_helper(x/2)
return (p*p+2*p*q,p*p+q*q)
Funciona porque puedes calcular fib(n),fib(n-1)
usando fib(n-1),fib(n-2)
si n es impar y si n es par, puedes calcular fib(n),fib(n-1)
usando fib(n/2),fib(n/2-1)
.
El caso base y el caso impar son simples. Para derivar el caso par, comience con a, b, c como valores consecutivos de fibonacci (por ejemplo, 8,5,3) y escríbalos en una matriz, con a = b + c. Darse cuenta:
[1 1] * [a b] = [a+b a]
[1 0] [b c] [a b]
A partir de eso, vemos que una matriz de los primeros tres números de fibonacci, multiplicada por una matriz de tres números consecutivos de fibonacci, es igual a la siguiente. Entonces sabemos que:
n
[1 1] = [fib(n+1) fib(n) ]
[1 0] [fib(n) fib(n-1)]
Asi que:
2n 2
[1 1] = [fib(n+1) fib(n) ]
[1 0] [fib(n) fib(n-1)]
Simplificar el lado derecho conduce al caso par.
¿Hay algún algoritmo para calcular el n-ésimo número de fibonacci en tiempo sub lineal?
Aquí hay una línea que calcula F (n), usando enteros de tamaño O (n), en O (log n) operaciones aritméticas:
for i in range(1, 50):
print(i, pow(2<<i, i, (4<<2*i)-(2<<i)-1)//(2<<i))
Usar números enteros de tamaño O (n) es razonable, ya que es comparable al tamaño de la respuesta.
Para entender esto, deje que phi sea la proporción áurea (la solución más grande para x ^ 2 = x + 1) y F (n) sea el n''th número de Fibonacci, donde F (0) = 0, F (1) = F (2) = 1
Ahora, phi ^ n = F (n-1) + F (n) phi.
Prueba por inducción: phi ^ 1 = 0 + 1 * phi = F (0) + F (1) phi. Y si phi ^ n = F (n-1) + F (n) phi, entonces phi ^ (n + 1) = F (n-1) phi + F (n) phi ^ 2 = F (n-1) phi + F (n) (phi + 1) = F (n) + (F (n) + F (n-1)) phi = F (n) + F (n + 1) phi. El único paso difícil en este cálculo es el que reemplaza a phi ^ 2 por (1 + phi), que sigue porque phi es la proporción áurea.
También los números de la forma (a + b * phi), donde a, b son números enteros se cierran bajo multiplicación.
Prueba: (p0 + p1 * phi) (q0 + q1 * phi) = p0q0 + (p0q1 + q1p0) phi + p1q1 * phi ^ 2 = p0q0 + (p0q1 + q1p0) phi + p1q1 * (phi + 1) = ( p0q0 + p1q1) + (p0q1 + q1p0 + p1q1) * phi.
Usando esta representación, se puede calcular phi ^ n en O (log n) operaciones enteras usando exponenciación por cuadratura. El resultado será F (n-1) + F (n) phi, desde el cual se puede leer el n''th número de Fibonacci.
def mul(p, q):
return p[0]*q[0]+p[1]*q[1], p[0]*q[1]+p[1]*q[0]+p[1]*q[1]
def pow(p, n):
r=1,0
while n:
if n&1: r=mul(r, p)
p=mul(p, p)
n=n>>1
return r
for i in range(1, 50):
print(i, pow((0, 1), i)[1])
Tenga en cuenta que la mayoría de este código es una función de exponenciación por cuadratura estándar.
Para llegar al trazador de líneas que inicia esta respuesta, se puede notar que representando phi por un número entero suficientemente grande X
, uno puede realizar (a+b*phi)(c+d*phi)
como la operación entera (a+bX)(c+dX) modulo (X^2-X-1)
. Entonces, la función de pow
puede ser reemplazada por la función de pow
Python estándar (que convenientemente incluye un tercer argumento z
que calcula el módulo de resultado z
. La X
elegida es 2<<i
.
De acuerdo con la referencia de Pillsy a la exponenciación de la matriz, tal que para la matriz
M = [1 1] [1 0]
entonces
fib(n) = Mn1,2
Elevar las matrices a los poderes usando la multiplicación repetida no es muy eficiente.
Dos enfoques para la exponenciación de la matriz son dividir y conquistar, lo que produce M n en los pasos O ( ln n ), o descomposición del valor propio que es tiempo constante, pero puede introducir errores debido a la precisión limitada del punto flotante.
Si desea un valor exacto mayor que la precisión de su implementación de coma flotante, debe usar el enfoque O (ln n) basado en esta relación:
Mn = (Mn/2)2 if n even = M·Mn-1 if n is odd
La descomposición de valores propios en M encuentra dos matrices U y Λ tales que Λ es diagonal y
M = U Λ U-1 Mn = ( U Λ U-1) n = U Λ U-1 U Λ U-1 U Λ U-1 ... n times = U Λ Λ Λ ... U-1 = U Λ n U-1 Elevar una matriz diagonal Λ a la n- ésima potencia es una cuestión simple de elevar cada elemento en Λ hasta la n- ésima, por lo que esto proporciona un método O (1) para elevar M a la n- ésima potencia. Sin embargo, los valores en Λ no son probablemente enteros, por lo que se producirá algún error.
Definiendo Λ para nuestra matriz 2x2 como
Λ = [ λ1 0 ] = [ 0 λ2 ]
Para encontrar cada λ , resolvemos
|M - λI| = 0
lo que da
|M - λI| = -λ ( 1 - λ ) - 1 λ² - λ - 1 = 0
usando la fórmula cuadrática
λ = ( -b ± √ ( b² - 4ac ) ) / 2a = ( 1 ± √5 ) / 2 { λ1, λ2 } = { Φ, 1-Φ } where Φ = ( 1 + √5 ) / 2
Si has leído la respuesta de Jason, puedes ver a dónde irá.
Resolviendo para los autovectores X 1 y X 2 :
if X1 = [ X1,1, X1,2 ] M.X1 1 = λ1X1 X1,1 + X1,2 = λ1 X1,1 X1,1 = λ1 X1,2 => X1 = [ Φ, 1 ] X2 = [ 1-Φ, 1 ]
Estos vectores dan U :
U = [ X1,1, X2,2 ] [ X1,1, X2,2 ] = [ Φ, 1-Φ ] [ 1, 1 ]
Invertir U usando
A = [ a b ] [ c d ] => A-1 = ( 1 / |A| ) [ d -b ] [ -c a ]
así que U -1 es dado por
U-1 = ( 1 / ( Φ - ( 1 - Φ ) ) [ 1 Φ-1 ] [ -1 Φ ] U-1 = ( √5 )-1 [ 1 Φ-1 ] [ -1 Φ ]
Prueba de cordura:
UΛU-1 = ( √5 )-1 [ Φ 1-Φ ] . [ Φ 0 ] . [ 1 Φ-1 ] [ 1 1 ] [ 0 1-Φ ] [ -1 Φ ] let Ψ = 1-Φ, the other eigenvalue as Φ is a root of λ²-λ-1=0 so -ΨΦ = Φ²-Φ = 1 and Ψ+Φ = 1 UΛU-1 = ( √5 )-1 [ Φ Ψ ] . [ Φ 0 ] . [ 1 -Ψ ] [ 1 1 ] [ 0 Ψ ] [ -1 Φ ] = ( √5 )-1 [ Φ Ψ ] . [ Φ -ΨΦ ] [ 1 1 ] [ -Ψ ΨΦ ] = ( √5 )-1 [ Φ Ψ ] . [ Φ 1 ] [ 1 1 ] [ -Ψ -1 ] = ( √5 )-1 [ Φ²-Ψ² Φ-Ψ ] [ Φ-Ψ 0 ] = [ Φ+Ψ 1 ] [ 1 0 ] = [ 1 1 ] [ 1 0 ] = M
Entonces el control de cordura se cumple.
Ahora tenemos todo lo que necesitamos para calcular M n 1,2 :
Mn = UΛnU-1 = ( √5 )-1 [ Φ Ψ ] . [ Φn 0 ] . [ 1 -Ψ ] [ 1 1 ] [ 0 Ψn ] [ -1 Φ ] = ( √5 )-1 [ Φ Ψ ] . [ Φn -ΨΦn ] [ 1 1 ] [ -Ψn ΨnΦ ] = ( √5 )-1 [ Φ Ψ ] . [ Φn Φn-1 ] [ 1 1 ] [ -Ψn -Ψn-1 ] as ΨΦ = -1 = ( √5 )-1 [ Φn+1-Ψn+1 Φn-Ψn ] [ Φn-Ψn Φn-1-Ψn-1 ]
asi que
fib(n) = Mn1,2 = ( Φn - (1-Φ)n ) / √5
Lo cual está de acuerdo con la fórmula dada en otra parte.
Puede derivarlo de una relación de recurrencia, pero en ingeniería computacional y simulación, el cálculo de valores propios y vectores propios de matrices grandes es una actividad importante, ya que proporciona estabilidad y armónicos de sistemas de ecuaciones, y permite elevar matrices a potencias altas de manera eficiente.
La aritmética de punto fijo es inexacta. El código C # de Jason da una respuesta incorrecta para n = 71 (308061521170130 en lugar de 308061521170129) y más.
Para la respuesta correcta, use un sistema de álgebra computacional. Sympy es una biblioteca para Python. Hay una consola interactiva en http://live.sympy.org/ . Copie y pegue esta función
phi = (1 + sqrt(5)) / 2
def f(n):
return floor(phi**n / sqrt(5) + 1/2)
Luego calcule
>>> f(10)
55
>>> f(71)
308061521170129
Puede intentar inspeccionar phi
.
Para los realmente grandes, esta función recursiva funciona. Utiliza las siguientes ecuaciones:
F(2n-1) = F(n-1)^2 + F(n)^2
F(2n) = (2*F(n-1) + F(n)) * F(n)
Necesitas una biblioteca que te permita trabajar con enteros grandes. Uso la biblioteca BigInteger desde https://mattmccutchen.net/bigint/ .
Comience con una serie de números de Fibonacci. Use fibs [0] = 0, fibs [1] = 1, fibs [2] = 1, fibs [3] = 2, fibs [4] = 3, etc. En este ejemplo, utilizo una matriz de los primeros 501 (contando 0). Puede encontrar los primeros 500 números de Fibonacci distintos de cero aquí: http://home.hiwaay.net/~jalison/Fib500.html . Se necesita una pequeña edición para ponerlo en el formato correcto, pero eso no es demasiado difícil.
Luego puede encontrar cualquier número de Fibonacci usando esta función (en C):
BigUnsigned GetFib(int numfib)
{
int n;
BigUnsigned x, y, fib;
if (numfib < 501) // Just get the Fibonacci number from the fibs array
{
fib=(stringToBigUnsigned(fibs[numfib]));
}
else if (numfib%2) // numfib is odd
{
n=(numfib+1)/2;
x=GetFib(n-1);
y=GetFib(n);
fib=((x*x)+(y*y));
}
else // numfib is even
{
n=numfib/2;
x=GetFib(n-1);
y=GetFib(n);
fib=(((big2*x)+y)*y);
}
return(fib);
}
He probado esto para el número 25,000 de Fibonacci y similares.
Puedes hacerlo exponiendo una matriz de enteros también. Si tienes la matriz
/ 1 1 /
M = | |
/ 1 0 /
entonces (M^n)[1, 2]
va a ser igual al n
° número de Fibonacci, si []
es un subíndice de la matriz y ^
es la exponenciación de la matriz. Para una matriz de tamaño fijo, la exponenciación a una potencia integral positiva se puede hacer en el tiempo O (log n) de la misma manera que con los números reales.
EDITAR: por supuesto, dependiendo del tipo de respuesta que desee, puede salirse con un algoritmo de tiempo constante. Al igual que las otras fórmulas muestran, el n
° número de Fibonacci crece exponencialmente con n
. Incluso con enteros sin signo de 64 bits, solo necesitará una tabla de búsqueda de 94 entradas para cubrir todo el rango.
SEGUNDA EDICIÓN: Hacer la matriz exponencial con una descomposición propia primero es exactamente equivalente a la solución de JDunkerly a continuación. Los valores propios de esta matriz son (1 + sqrt(5))/2
y (1 - sqrt(5))/2
.
Puedes usar la extraña ecuación de raíz cuadrada para obtener una respuesta exacta. La razón es que $ / sqrt (5) $ cae al final, solo tienes que hacer un seguimiento de los coeficientes con tu propio formato de multiplicación.
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
Uno de los ejercicios en SICP es sobre esto, que tiene la respuesta que se describe here.
En el estilo imperativo, el programa se vería como
Function Fib(count) a ← 1 b ← 0 p ← 0 q ← 1 While count > 0 Do If Even(count) Then p ← p² + q² q ← 2pq + q² count ← count ÷ 2 Else a ← bq + aq + ap b ← bp + aq count ← count - 1 End If End While Return b End Function
Wikipedia tiene una solución de formulario cerrado http://en.wikipedia.org/wiki/Fibonacci_number
O en c #:
public static int Fibonacci(int N)
{
double sqrt5 = Math.Sqrt(5);
double phi = (1 + sqrt5) / 2.0;
double fn = (Math.Pow(phi, N) - Math.Pow(1 - phi, N)) / sqrt5;
return (int)fn;
}
usando R
l1 <- (1+sqrt(5))/2
l2 <- (1-sqrt(5))/2
P <- matrix(c(0,1,1,0),nrow=2) #permutation matrix
S <- matrix(c(l1,1,l2,1),nrow=2)
L <- matrix(c(l1,0,0,l2),nrow=2)
C <- c(-1/(l2-l1),1/(l2-l1))
k<-20 ; (S %*% L^k %*% C)[2]
[1] 6765
ver el algoritmo de divide y vencerás here
El enlace tiene pseudocódigo para la exponenciación de la matriz mencionada en algunas de las otras respuestas para esta pregunta.
Si quieres el número exacto (que es un "bignum", en lugar de un int / float), entonces me temo que
¡Es imposible!
Como se indicó anteriormente, la fórmula para los números de Fibonacci es:
fib n = piso (phi n / √5 + 1/2)
fib n ~ = phi n / √5
¿Cuántos dígitos es fib n
?
numDigits (fib n) = log (fib n) = log (phi n / √5) = log phi n - log √5 = n * log phi - log √5
numDigits (fib n) = n * const + const
es O ( n )
Como el resultado solicitado es de O ( n ), no se puede calcular en menos de O ( n ) tiempo.
Si solo desea los dígitos más bajos de la respuesta, entonces es posible calcular en tiempo sub-lineal usando el método de exponenciación de la matriz.
El n
° número de Fibonacci está dado por
f(n) = Floor(phi^n / sqrt(5) + 1/2)
dónde
phi = (1 + sqrt(5)) / 2
Suponiendo que las operaciones matemáticas primitivas ( +
, -
, *
y /
) son O(1)
puede usar este resultado para calcular el n
° número de Fibonacci en el tiempo O(log n)
( O(log n)
debido a la exponenciación en la formula).
Cª#:
static double inverseSqrt5 = 1 / Math.Sqrt(5);
static double phi = (1 + Math.Sqrt(5)) / 2;
/* should use
const double inverseSqrt5 = 0.44721359549995793928183473374626
const double phi = 1.6180339887498948482045868343656
*/
static int Fibonacci(int n) {
return (int)Math.Floor(Math.Pow(phi, n) * inverseSqrt5 + 0.5);
}