Algoritmos para calcular la sucesión de Fibonnaci usando latexify para generar el código LaTeX

En este artículo uso como excusa la sucesión de Fibonacci y algunos algoritmos que permiten su cálculo para analizar la forma en que puedo usar el paquete latexify de python para obtener, a partir del código python y de forma directa, los algoritmos de las funciones así como las expresiones de las funciones matemáticas en código LaTeX. Además uso python para realizar los cálculos, pandas para realizar las tablas y el paquete pythontex para crear de forma dinámica el pdf resultado de compilar el fichero LyX. Los algoritmos usados se pueden encontrar por internet en las páginas referenciadas. Termino el artículo con el cálculo del término 10000 de la sucesión.

La sucesión de Fibonacci es la sucesión infinita de números naturales

donde el primer elemento es 0, el segundo es 1 y cada elemento restante es la suma de los dos anteriores. A cada elemento de esta sucesión se le llama número de Fibonacci. Se debe a Leonardo de Pisa (1170-1250), matemático italiano del siglo XIII también conocido como Fibonacci (hijo de Bonacci [1]). La ideó como la solución a un problema de la cría de conejos:

Cierto hombre tenía una pareja de conejos juntos en un lugar cerrado y uno desea saber cuántos son creados a partir de este par en un año cuando es su naturaleza parir otro par en un simple mes, y en el segundo mes los nacidos parir también.

Para obtener sus términos, se cumple la relación de recurrencia siguiente:

$\mathrm{fib}(x)=\left\{ \begin{array}{ll} 0,$& \mathrm{if}\ x=0\\ 1, & \mathrm{if}\ x=1\\ \mathrm{fib}\mathopen{}\left(x-1\mathclose{}\right)+\mathrm{fib}\mathopen{}\left(x-2\mathclose{}\right), & \mathrm{otherwise} \end{array}\right.

Podemos programar fácilmente el algoritmo que se deduce de la fórmula anterior a partir de una función recursiva:

Enpython quedaría:

def fib(x):
  if x == 0:
    return 0
  elif x == 1:
    return 1
  else:
    return fib(x-1) + fib(x-2)

con la función anterior, por ejemplo, podemos obtener que $fib(20)=6765$. Pero su eficiencia es muy pobre. La programación recursiva de la función de Fibonacci tiene una complejidad, como mínimo exponencial (véase, por ejemplo https://www.ugr.es/~eaznar/fibo.htm). De hecho, calcular el término 100 con la función anterior ya no es factible.

Existen algoritmos mucho más eficientes (por ejemplo la solución iterativa aunque la forma directa tampoco es buena) y me voy a referir a dos de ellos, el primero tomado de la Wikipedia https://es.wikipedia.org/wiki/Sucesi%C3%B3n_de_Fibonacci. Para construirlo debemos usar dos cuestiones:

1. Usar un algoritmo óptimo para calcular potencias:

   $$x^{n}=\begin{cases} 1$$& n=0\\ \left(x^{n/2}\right)^{2} & n\,\,{\textstyle par}\\ x\cdot x^{n-1} & n\,\,{\textstyle impar} \end{cases}

   o equivalentemente:

   $\mathrm{potencia}(x,n)=\left\{ \begin{array}{ll} x,$& \mathrm{if}\ n=1\\ \mathopen{}\left(\mathrm{potencia}\mathopen{}\left(x,\left\lfloor \frac{n}{2}\right\rfloor \mathclose{}\right)\mathclose{}\right)^{2}, & \mathrm{if}\ n\mathbin{\%}2=0\\ x\cdot\mathrm{potencia}\mathopen{}\left(x,n-1\mathclose{}\right), & \mathrm{otherwise} \end{array}\right.

   cuyo algoritmo es:

   y enpython podría ser:

   # Potencia
   def potencia(x,n):
       if n==1:
           return x
       elif n % 2 == 0:
           return (potencia(x,n // 2))**2
       else:
           return x*potencia(x,n-1)
   

   Pero no vamos a usar el algoritmo anterior, sino otro (véase, por ejemplo, https://www.geeksforgeeks.org/fast-exponentiation-in-python/) para calcular la potencia de un número de forma iterativa que consiste en

   y enpython podría ser:

   # Potencia
   def pot_ite(x,n):
       i=n
       a = 1
       c = x
       while i > 0:
           if i % 2 != 0:
               a = a * c
           c = c**2
           i = i // 2
       return a
   

   Con él podemos obtener, por ejemplo, que pot_ite(2,10) $=2^{10}=1024$

2. Además, podemos usar matrices para determinar los términos de la sucesión de Fibonacci

   Proposición

   $$\left[\begin{array}{cc} 1$$& 1\\ 1 & 0 \end{array}\right]^{n}=\left[\begin{array}{cc} F(n+1) & F(n)\\ F(n) & F(n-1) \end{array}\right]

   Demostración

   Se demuestra por inducción:

   • $n=1$

     $\left[\begin{array}{cc} 1$& 1\\ 1 & 0 \end{array}\right]^{1}=\left[\begin{array}{cc} F(2) & F(1)\\ F(1) & F(0) \end{array}\right] que es cierto

   • Suponemos cierto para $n\longmapsto$ $\left[\begin{array}{cc} 1$& 1\\ 1 & 0 \end{array}\right]^{n}=\left[\begin{array}{cc} F(n+1) & F(n)\\ F(n) & F(n-1) \end{array}\right]

   • Veamos que es cierto para $n+1$

     $$\begin{aligned} \left[\begin{array}{cc} 1$$& 1\\ 1 & 0 \end{array}\right]^{n+1} & =\left[\begin{array}{cc} 1 & 1\\ 1 & 0 \end{array}\right]^{n}\left[\begin{array}{cc} 1 & 1\\ 1 & 0 \end{array}\right]=\\&=\left[\begin{array}{cc} F(n+1) & F(n)\\ F(n) & F(n-1) \end{array}\right]\left[\begin{array}{cc} 1 & 1\\ 1 & 0 \end{array}\right]=\\ = & \left[\begin{array}{cc} F(n+1)+F(n) & F(n+1)+0\\ F(n)+F(n-1) & F(n)+0 \end{array}\right]=\\&=\left[\begin{array}{cc} F(n+2) & F(n+1)\\ F(n+1) & F(n) \end{array}\right]\end{aligned}

Observar que las matrices

quedan completamente determinada solo por dos valores, si $F(n-1)=a$ y $F(n)=b\Rightarrow F(n+1)=a+b$, es decir, todas son de la forma $\left[\begin{matrix}a+b$& b\\ b & a \end{matrix}\right]

Si llamamos $c=F(k-1)$ y $d=F(k)\,\Rightarrow\,F(k+1)=c+d$, por tanto

Es decir,

• si $n=2\cdot k$ ($n$ par) y conozco $F(k-1)=c$ y $F(k)=d$ puedo calcular, a partir de los valores anteriores, tanto $F(n-1)=F(2k-1)=c^{2}+d^{2}$ como $F(n)=F(2k)=d\cdot(2\cdot c+d)$

• si $n=2\cdot k+1$ ($n$ es impar) y aplico el algoritmo Potencia iterativa tendría que multiplicar las matrices

  $$A*C=\left[\begin{matrix}a+b$$& b\\ b & a \end{matrix}\right]\cdot\left[\begin{matrix}c+d & d\\ d & c \end{matrix}\right]=\left[\begin{matrix}bd+\left(a+b\right)\left(c+d\right) & bc+d\left(a+b\right)\\ ad+b\left(c+d\right) & ac+bd \end{matrix}\right]\space{}

  y en consecuencia $F(n)=b\cdot c+d\cdot(a+b)$ y $F(n-1)=a\cdot c+b\cdot d$

A partir de lo anterior y el algoritmo para calcular la potencia tenemos que:

y en python:

# Solución de wikipedia
def fib_optima(n):
    if n<=0:
        return 0
    # Llegamos hasta el n-1 ya que el último (F(n)) se obtiene sumando a+b
    i=n-1
    a, b = 1, 0
    c, d = 0, 1
    while i > 0:
        if i % 2 != 0:
            a, b = a*c+b*d, b*c+d*(a+b)
        c, d = c**2+d**2, d*(2*c+d)
        i = i // 2
    return a+b

con él, por ejemplo, ya sí podemos calcular:

o:

En Project Nayuki tenemos otra forma de calcular los términos de la sucesión de forma incluso más eficiente, en este caso nos basamos en que:

Y en consecuencia obtenemos que:

Es decir, si conozco $F(k)=a$ y $F(k+1)=b$ puedo calcular $c=F(2k)=a\cdot(2\cdot b-a)$ y $d=F(2k+1)=a^{2}+b^{2}$, por tanto:

• Si $n=2\cdot k$ es par $\Rightarrow F(n)=F(2\cdot k)=c=a\cdot(2\cdot b-a)$ y $F(n+1)=F(2k+1)=d=a^{2}+b^{2}$
• Si $n=2\cdot k+1$ es impar $F(n)=F(2\cdot k+1)=d$ y $F(n+1)=F(2k+2)=F(2k)+F(2k+1)=c+d$

y obtenemos los algoritmos

y sus correspondientes funciones en python:

# Algoritmo de https://www.nayuki.io/page/fast-fibonacci-algorithms
# (Pública) Retorna F(n).
def fibonacci(n):
    if n < 0:
        return 0
    return _fibo(n)[0]


# (Privada) Retorna la tupla (F(n), F(n+1)).
def _fibo(n):
    if n == 0:
        return (0, 1)
    else:
        a, b = _fibo(n // 2)
        c = a * (2 * b - a)
        d = a**2 + b**2
        if n % 2 == 0:
            return (c, d)
        else:
            return (d, c + d)

A partir de esta última, por ejemplo, podemos plantear y comprobar las cuestiones:

1. Halla el cociente $\frac{\textrm{fi}(n+1)}{\textrm{fib}(n)}$ para los $45$ primeros términos de la sucesión (salvo el primero). ¿A qué número famoso se parece cada vez más?

   n	F(n)	F(n+1)	Cociente
   1	1	1	1.0000000000000000
   2	1	2	2.0000000000000000
   3	2	3	1.5000000000000000
   4	3	5	1.6666666666666667
   5	5	8	1.6000000000000001
   6	8	13	1.6250000000000000
   7	13	21	1.6153846153846154
   8	21	34	1.6190476190476191
   9	34	55	1.6176470588235294
   10	55	89	1.6181818181818182
   11	89	144	1.6179775280898876
   12	144	233	1.6180555555555556
   13	233	377	1.6180257510729614
   14	377	610	1.6180371352785146
   15	610	987	1.6180327868852460
   16	987	1597	1.6180344478216819
   17	1597	2584	1.6180338134001253
   18	2584	4181	1.6180340557275541
   19	4181	6765	1.6180339631667064
   20	6765	10946	1.6180339985218033
   21	10946	17711	1.6180339850173580
   22	17711	28657	1.6180339901755971
   23	28657	46368	1.6180339882053250
   24	46368	75025	1.6180339889579021
   25	75025	121393	1.6180339886704431
   26	121393	196418	1.6180339887802426
   27	196418	317811	1.6180339887383031
   28	317811	514229	1.6180339887543225
   29	514229	832040	1.6180339887482036
   30	832040	1346269	1.6180339887505408
   31	1346269	2178309	1.6180339887496482
   32	2178309	3524578	1.6180339887499890
   33	3524578	5702887	1.6180339887498589
   34	5702887	9227465	1.6180339887499087
   35	9227465	14930352	1.6180339887498896
   36	14930352	24157817	1.6180339887498969
   37	24157817	39088169	1.6180339887498940
   38	39088169	63245986	1.6180339887498951
   39	63245986	102334155	1.6180339887498947
   40	102334155	165580141	1.6180339887498949
   41	165580141	267914296	1.6180339887498949
   42	267914296	433494437	1.6180339887498949
   43	433494437	701408733	1.6180339887498949
   44	701408733	1134903170	1.6180339887498949
   45	1134903170	1836311903	1.6180339887498949

   es decir: $\lim_{n\to\infty}\frac{\textrm{fib}(n+1)}{\textrm{fib}(n)}=\varphi=\dfrac{1+\sqrt{5}}{2}\approx1.6180339887498948482045\ldots$

2. Halla el máximo común divisor de dos números de Fibonacci ¿qué se obtiene? ¿existe alguna relación entre ellos?

   n	m	F(n)	F(m)	MCD(n,m)	F(MCD(n,m))	MCD(F(n),F(m))
   10	20	55	6765	10	55	55
   15	31	610	1346269	1	1	1
   20	31	6765	1346269	1	1	1
   25	52	75025	32951280099	1	1	1
   30	39	832040	63245986	3	2	2
   35	52	9227465	32951280099	1	1	1
   40	53	102334155	53316291173	1	1	1
   45	56	1134903170	225851433717	1	1	1
   50	69	12586269025	117669030460994	1	1	1
   55	63	139583862445	6557470319842	1	1	1
   60	88	1548008755920	1100087778366101931	4	3	3
   65	88	17167680177565	1100087778366101931	1	1	1
   70	90	190392490709135	2880067194370816120	10	55	55
   75	84	2111485077978050	160500643816367088	3	2	2

   es decir, el máximo común divisor de dos números de Fibonacci es otro número de Fibonacci y además $MCD(F(n),F(m))=F(MCD(n,m))$MCD(F(n),F(m))=F(MCD(n,m)).

   De lo anterior se deduce que ${\displaystyle F(n)}$ y $F(n+1)$ son primos relativos y que $F(k)$ divide exactamente a $F(n\cdot k)$

Para terminar, calculemos con el último algoritmo el término $10000$ de la sucesión de Fibonacci:

con un tiempo de ejecución para calcularlo de 2.9325485229492188e-05 segundos.

Fichero fuente y el pdf final de una posible compilación.

• Fichero en formato LyX con el código python necesario para hacerlo.
• Fichero en formato pdf que se obtiene como resultado de compilar el fichero anterior.