Números de Fibonacci, su complejidad y su programación.

La llamada sucesión de los números de Fibonacci
1, 1, 2, 3, 5, 8, 13, 21, 34, 55, ...
que aparece por primera vez en el célebre, Liber abaci (1202), publicado por Leonardo Pisano, como solución a un problema de reproducción de conejos, satisface la ecuación en recurrencia lineal
F_n = F_n-1 + F_n-2
donde F₁ = 1 = F₂

Cláramente, no es la única sucesión de números, solución de la misma ecuación en recurrencia. Por ejemplo, la sucesión
1, 3, 4, 7, 11, 18, 29, 47, 76, ...
llamados los números de Lucas (matemático francés, 1842-1891). También satisfacen la misma ecuación en recurrencia. Pero la diferencia está en los dos primeros términos, que en este caso son F₁ = 1 ; F₂ = 3.

Como los siguientes valores de una tal sucesión F₁, F₂, F₃, ... están determinados por los dos primeros. Hay tantas sucesiones solución como parejas de números F₁ = a ; F₂ = b.

Programación recursiva

Si queremos programar una función que calcule el n-ésimo término, F_n, de una de estas soluciones. La primera aproximación es precísamente recursiva. Por ejemplo, hoy día, casi todos los lenguajes de alto nivel, admiten las tres siguientes líneas que implementan el cálculo de los números de Fibonacci, con pocas diferencias de sintaxis

F(1) = 1;

F(2) = 1;

F(n) = F(n-1) + F(n-2);

El problema es su complejidad. No sólo en espacio. Ya que todas las llamadas anteriores, a la propia función, han de recordarse, calcularse y asignarse para poder hallar F(n). Sino también en número de operaciones elementales. Que se traduce en complejidad de tiempo. Ya que se suman todas las operaciones elementales necesarias para calcular todos los anteriores.

Veamos una aproximación ingenua, pero precisa, al cálculo de la complejidad de la función recursiva anterior. Si llamamos S(n) al número de sumas necesarias para hallar F(n). Para los primeros valores se tiene
S(1) = 0 = S(2), S(3)=1, S(4) = 2, S(5) = 4, S(6) = 7, ...
y en general por inducción, el número de sumas para calcular F(n) es igual a
S(n) = S(n-1) + S(n-2) + 1
Su crecimiento es mas rápido que el la función original. Ya que se obtiene cláramente por inducción que
F(n-2) < S(n)

Complejidad.

Pero ¿Cómo de rápido crece la función de Fibonacci?.
Para responder a esa pregunta. Calcularemos las raíces de la ecuación característica asociada a la ecuación en recurrencia F(n) = F(n-1) + F(n-2).
x² = x + 1
o sea, las raíces de x² - x - 1 = 0

Una de ellas, es el llamado número de oro (Golden Ratio), cuyo valor aproximado es 1.61803 y su valor exacto es c = (1+√5)/2.

Cláramente, se verifica c² = c + 1 y por tanto, para todo n mayor que dos, también
cⁿ = c^n-1 + c^n-2

O sea, la progresión geométrica {cⁿ} satisface la misma ecuación en recurrencia que la función de Fibonacci F(n) = F(n-1) + F(n-2). Ahora, por inducción

Como, c⁰ = 1 = F(2), c = c¹ < 2 = F(3).
obtenemos que c^n-2 < F(n)

En conclusión, la función de Fibonacci crece, como mínimo, exponencialmente. Ahora, enlazando las dos desigualdades, para todo n se tiene

c^n-4 < F(n-2) < S(n)

y también la función S(n) crece, como mínimo, exponencialmente.

Por tanto, la programación recursiva de la función de Fibonacci tiene una complejidad, como mínimo, exponencial. Y eso, independientemente, de lo bien que gestiene el compilador o intérprete correspondiente la programación recursiva. Veamos a continuación que una programación iterativa, sin embargo, tiene una complejidad mucho mejor.

Programación iterativa

Para calcular iterativamente la sucesión de Fibonacci, F_n, basta con las siguiente líneas de asignación, que con pocas diferencias de sintaxis, admiten casi todos los lenguajes de alto nivel

lista (a,b) = (1,1);

for(i=2;i<=n;i++){(a,b) = (b, a+b);}

return b;

Cláramente, se asignan sólo dos variables en cada iteración y se hace también sólo una suma. Por tanto, el número de sumas totales ahora es
S(n) = n
y su complejidad operacional, en este caso, es lineal, O(n), en vez de exponencial.

Mejorando la programación iterativa

Si se considera el producto de matrices

Es fácil desmostrar por inducción que las potencias de la matriz están relacionadas doblemente con la sucesión de los números de Fibonacci, por la igualdad

Aplicando por tanto el algoritmo de la exponenciación rápida a la matriz es posible calcular el término n-ésimo, F_n de la sucesión de Fibonacci, en log₂(n) iteraciones.

Una primera implementación en el paquete Mathematica de Wolfram Research sería

          Fibo3[n_Integer?Positive] := 
	             Module[{ mFibo={{1,1},{1,0}}, mRes={{1,0},{0,1}}, contador=n-1}, 
	               While[contador>1,
                         If[OddQ[contador], mRes=mRes.mFibo];
                         contador=Floor[contador/2];
                         mFibo=mFibo.mFibo
                       ];
                       mFibo[[1,1]] mRes[[1,1]] + mFibo[[2,1]] mRes[[1,2]]
                     ]

Como en cada iteración se hacen como máximo 2 multiplicaciones de matrices 2x2. Existe una constante λ tal que la función de complejidad operacional estará acotada, S(n) ≤ λlog₂(n).

Obteniendose una complejidad logarítmica, O(log₂(n)) para esta versión del algoritmo.

Además, como la complejidad de las multiplicación es mayor que la la suma, podemos medir la complejidad operacional contando el número de multiplicaciones elementales. Como para calcular el producto de dos matrices 2x2, hace falta 8 multiplicaciones de sus entradas. Y en cada iteración puede haber dos multiplicaciones de matrices, como mucho se efectúan 16 multiplicaciones de elementos (números enteros positivos). De donde λ = 16.

Mejorando todavía más.

Todavía se puede mejorar la eficacia de la implementación, teniendo en cuenta que las potencias sucesivas de la matriz anterior, o sea todas las matrices implicadas, salen todas simétricas y bastarían 3 elementos para determinarla. Además, como también todas verifican la relación invariante a₁₁ - a₂₂ = a₁₂. En realidad, bastan sólo dos elementos para determinar el producto.

Una programación que tenga en cuenta esto, mejoraría la constante hasta la mitad, λ = 8, de la función complejidad última. Aunque seguiría siendo logarítmica. Por ejemplo en Mathematica.

          Fibo4[n_Integer?Positive]:=
	             Module[{m11=1, m12=1, m22=0, r11=1, r12=0, r22=1,contador=n-1},
	                While[contador>1,
                          If[OddQ[contador],{r11, r22}= m12*r12 + {m11*r11, m22*r22}; 
                          r12 = r11 - r22;];
                          contador=Floor[contador/2];
                          {m11, m22} = m12^2 + {m11^2, m22^2}; m12 = m11 - m22;
                        ]; 
                        m11 r11 + m12 r12
                      ]

Comparación con la fórmula de Binet.

Desde 1841, sabemos que el término n-ésimo de la sucesión de Fibonacci, puede calcularse por la llamada fórmula de Binet (Jacques Binet, 1786 - 1856). Aunque, en realidad el primero en demostrarla fue Abraham de Moivre (1667 - 1754).

La programación exacta, o sea sin calcular ninguna raíz cuadrada, sólo desarrollando formalmente las potencias y cancelando hasta conseguir un número natural, es posible hacerla de forma sencilla en algún lenguaje que tenga implentado dicho proceso. Por ejemplo en Mathematica se puede hacer con el siguiente código:

          FiboBinet[n_Integer?Positive]:=(((1 + Sqrt[5])/2)^n - ((1 - Sqrt[5])/2)^n)/Sqrt[5] // Expand

En este caso, aunque la exponenciación rápida se aplica en el desarrollo de las dos potencias, obteniendose un coste igual que el último anterior, hay que sumar el coste de la función Expand (cancelar hasta dar una forma normal) que lo hacen mas costoso y mas lento. Una prueba, en Mathematica con Timing, demuestra claramente la superioridad de la anterior implemantación con respecto a esta.

Si por el contrario, se calcula un valor previo de los irracionales implicados √5, (1 + √5)/2, (1 - √5)/2 con un número de decimales fijos, la fórmula de binet tiene la misma complejidad que el mejor algoritmo iterativo, 8log₂(n), pero no sirve para todo n. Ya que da resultados inexactos a partir de un n, que depende de la precisión inicial, en adelante.

Una implementación en Mathematica, que da resultados correctos hasta n = 19999, con output números con más de 4.000 cifras decimales, puede ser

             a=N[5^(-1/2),10^4];
             b=N[GoldenRatio,10^4];
             moivre[n_Integer?(#>=0&)]:=If[Mod[n,2]==0,IntegerPart[b^n*a],IntegerPart[b^n*a+1]]

Epílogo y conclusiones.

Finalmente, observamos que la programación o implementación recursiva dista mucho de ser eficaz. Mientras que la fórmula explícita es muy elegante pero inexacta en la práctica a partir de un n en adelante. E incluso ineficaz cuando se compara con un método iterativo específicamente diseñado.