Sucesion de Fibonacci recursiva optimizada

[do-while]

¡Buenas!

Esta mañana me he dedicado a perder el tiempo, así que comparto el tiempo perdido con vosotros y no tengo la sensación de haberlo perdido para nada.

Se trata de la sucesión de Fibonacci (ya conocida por todos (supongo)).

El algoritmo obtenido no es tan rápido como la versión iterativa (creo que anda cerca), pero mejora con creces la versión recursiva.

Evidentemente no es inmediata, como la formula que obtiene el enésimo termino en función de n y de phi (el numero aureo...)

Tampoco he utilizado programación dinamica, así que el único problema de memoria, si lo hubiese, estaría originado por una cantidad excesiva de llamadas recursivas.

Aquí os dejo un par de códigos. El primero con explicaciones sobre como se llega a la formula utilizada y con una comparación entre el algoritmo clásico y el optimizado, y el segundo para utilizarlo desde la linea de comandos, por si a alguien le hace ilusión.

¡Saludos!

Código (c) [Seleccionar] Expandir


/*
  JUGANDO UN POCO CON LA SUCESION DE FIBONACCI
  ============================================

  La sucesion de Fibonacci se define de forma recursiva como:

fibo(n) = n  si n = 0 ó 1
fibo(n) = fibo(n - 1) + fibo(n - 2) si n > 1

  Aplicar la recursion de esta forma hace que el numero de llamadas sea de orden
exponencial con respecto al tamaño inicial del dato, por lo que el algoritmo se
vuelve lento rapidamente.

  Suponiendo n suficientemente grande podemos descomponer fibo(n - 1) en el paso
de recurcion como sigue:

fibo(n) = fibo(n - 1) + fibo(n - 2) = (fibo(n - 2) + fibo(n - 3)) + fibo(n - 2) =
= 2 * fibo(n - 2) + fibo(n - 3)

  Aplicamos el mismo proceso a fibo(n - 2):

fibo(n) = 2 * fibo(n - 2) + fibo(n - 3) = 2 * (fibo(n - 3) + fibo(n - 4)) + fibo(n - 3) =
= 3 * fibo(n - 3) + 2 * fibo(n - 4)

  Y para ver las cosas mas claras lo repetimos con fibo(n - 3):

fibo(n) = 3 * fibo(n - 3) + 2 * fibo(n - 4) = 3 * (fibo(n - 4) + fibo(n - 5) + 2 * fibo(n - 4) =
= 5 * fibo(n - 4) + 3 * fibo(n - 5)

  Por si no queda claro donde queremos llegar resumimos lo que tenemos hasta ahora:

  Cada una de las siguientes expresion son iguales a fibo(n):

1 * fibo(n - 1) + 1 * fibo(n - 2) = fibo(2) * fibo(n - 1) + fibo(1) * fibo(n - 2) =
2 * fibo(n - 2) + 1 * fibo(n - 3) = fibo(3) * fibo(n - 2) + fibo(2) * fibo(n - 3) =
3 * fibo(n - 3) + 2 * fibo(n - 4) = fibo(4) * fibo(n - 3) + fibo(3) * fibo(n - 4) =
5 * fibo(n - 4) + 3 * fibo(n - 5) = fibo(5) * fibo(n - 4) + fibo(4) * fibo(n - 5)

  Si seguimos de esta forma podemos intuir que

  fibo(n) = fibo(k + 1) * fibo(n - k) + fibo(k) * fibo(n - k - 1) := F(n,k) para 0 <= k < n

  Se puede demostrar por induccion que fijando un N natural (N >= 2) se tiene que

  fibo(N) = F(N,k) para todo 0 <= k < N

  y que fijado K natural para todo n > K se tiene que

  fibo(n) = F(n,K)

  por lo tanto tenemos:

  fibo(n) = F(n,k) = fibo(k + 1) * fibo(n - k) + fibo(k) * fibo(n - k - 1) para 0 <= k < n


  SUCESION DE FIBONACCI OPTIMIZADA
  ================================

  Escribiremos fibo(n) como  f(n), es mas corto :-)
  Partiendo de la expresion F(n,k) consideraremos dos casos:

  - n par:

        Sea k = n / 2

        Asi F(n,k) = f(k + 1) f(n - k) + f(k) f(n - k - 1) =
        f(n/2 + 1) f(n/2) + f(n/2) f(n/2 - 1) =
        =f(n/2) (f(n/2 + 1) + f(n/2 - 1)

        F(n,n/2) = f(n / 2) * (f(n / 2 + 1) + f(n / 2 - 1)) := F(n)

    - n impar:

        Sea k = (n - 1) / 2

        F(n,k) = f(k + 1) f(n - k) + f(k) f(n - k - 1) =
        f((n - 1) / 2 + 1) f((n + 1) / 2) + f((n - 1) / 2) f((n + 1) / 2 - 1) =
        f((n + 1) / 2) f((n + 1) / 2) + f((n - 1) / 2) f((n - 1) / 2) =
        (f((n + 1) / 2))^2 + (f((n - 1) / 2))^2

        F(n,(n - 1)/2) = (f((n + 1) / 2))^2 + (f((n - 1) / 2))^2 := F(n)

    Por lo tanto definimos la sucesion de Fibonacci como sigue:

    fibo(0) = 0
    fibo(1) = 1
    fibo(n) = F(n) si n > 1

    donde:
        F(n) = fibo(n / 2) * (fibo(n / 2 + 1) + fibo(n / 2 - 1)) si n par
        F(n) = (fibo((n + 1) / 2))^2 + (fibo((n - 1) / 2))^2     si n impar

    Observemos ahora que si n==2,

    fibo(2) = F(2) = fibo(1) * (fibo(2) + fibo(0)) = fibo(1) * fibo(2)

    y tenemos una recursion infinita, ya que para calcular fibo(2) hay que calcular fibo(2)...

    Asi pues, tenemos la definicion final de la sucesion de Fibonacci:

    fibo(0) = 0
    fibo(1) = 1
    fibo(2) = 1
    fibo(n) = F(n) si n > 1

    con F(n) como antes.

    Definiendo fibo(n) de esta forma, en cada paso de recursion generamos dos o tres llamadas
    a la misma funcion fibo, pero en cada una de ellas reducimos el tamaño del dato inicial a
    la mitad, por lo que el coste computacional es mucho menor (¿Cual es el coste? No se calcularlo)
*/

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define MIN 37
#define MAX 42

typedef unsigned long long fibo_t;

unsigned long long llamadas_clasico;
unsigned long long llamadas_optimizado;

fibo_t fibo_clasico(fibo_t n)
{
    llamadas_clasico++;

    if(n == 0)
        return 0;
    if(n == 1)
        return 1;

    return fibo_clasico(n - 1) + fibo_clasico(n - 2);
}

fibo_t fibo(fibo_t n)
{
    llamadas_optimizado++;

    if(n == 0)
        return 0;
    if(n == 1)
        return 1;
    if(n == 2) //Evitamos recurion infinita si n == 2
        return 1;

    if(n % 2) //n impar y mayor que 2
    {
        fibo_t k1,k2;

        k1 = fibo((n + 1) / 2);
        k2 = fibo((n - 1) / 2);

        return k1 * k1 + k2 * k2;
    }

    //n par y mayor que 2
    return fibo(n / 2) * (fibo((n / 2) + 1) + fibo((n / 2) - 1));
    //si n == 2 fibo(n / 2) * fibo(n / 2 + 1) = fibo(1) * fibo(2)
    //habria recursion infinita si no se pone n == 2 como caso base
}

int main(int argc, char *argv[])
{
    fibo_t i,valor;
    time_t ini;

    printf("Fibonacci clasico:\n");

    ini = time(NULL);

    for(i = MIN ; i < MAX ; i++)
    {
        llamadas_clasico = 0;
        valor = fibo_clasico(i);

        printf("  fibo_clasico(%llu) = %llu (%llu llamadas)\n",i,valor,llamadas_clasico);
    }

    printf("%llu segundos\n",time(NULL) - ini);

    printf("\nFibonacci optimizado:\n");

    ini = time(NULL);

    for(i = MIN ; i < MAX ; i++)
    {
        llamadas_optimizado = 0;
        valor = fibo(i);

        printf("  fibo(%llu) = %llu (%llu llamadas)\n",i,valor,llamadas_optimizado);
    }

    printf("%llu segundos\n",time(NULL) - ini);

    return 0;
}


Código (c) [Seleccionar] Expandir


#include <stdio.h>

typedef unsigned long long fibo_t;

fibo_t fibo(fibo_t n)
{
    if(n == 0)
        return 0;
    if(n == 1)
        return 1;
    if(n == 2)
        return 1;

    if(n % 2)
    {
        fibo_t k1,k2;

        k1 = fibo((n + 1) / 2);
        k2 = fibo((n - 1) / 2);

        return k1 * k1 + k2 * k2;
    }

    return fibo(n / 2) * (fibo((n / 2) + 1) + fibo((n / 2) - 1));
}

int convertir(char *s, fibo_t *valor)
{
    *valor = 0;

    while(*s)
    {
        if((*s) < '0' || (*s) > '9')
        {
            *valor = 0;
            return 0;
        }

        *valor *= 10;
        *valor += (*(s++)) - '0';
    }

    return 1;
}

int main(int argc, char *argv[])
{
    int i;
    fibo_t n;

    if(argc < 2)
    {
        fprintf(stderr,"\nLinea de comandos: %s PARAMETRO_1 [PARAMETRO_2 ... PARAMETRO_N]\n\n",argv[0]);
        fprintf(stderr,"Los distintos parametros deben ser enteros positivos\n");
        return 1;
    }

    fprintf(stdout,"\n");

    for(i = 1 ; i < argc ; i++)
    {
        if(!convertir(argv[i],&n))
            fprintf(stderr,"%s: parametro incorrecto\n",argv[i]);
        else
            fprintf(stdout,"%s(%llu) = %llu\n",argv[0],n,fibo(n));
    }

    return 0;
}


[0xDani]

Muy buen código, y de hecho si que está optimizado 

Fibonacci clasico:
  fibo_clasico(37) = 24157817 (78176337 llamadas)
  fibo_clasico(38) = 39088169 (126491971 llamadas)
  fibo_clasico(39) = 63245986 (204668309 llamadas)
  fibo_clasico(40) = 102334155 (331160281 llamadas)
  fibo_clasico(41) = 165580141 (535828591 llamadas)
46 segundos

Fibonacci optimizado:
  fibo(37) = 24157817 (88 llamadas)
  fibo(38) = 39088169 (141 llamadas)
  fibo(39) = 63245986 (90 llamadas)
  fibo(40) = 102334155 (131 llamadas)
  fibo(41) = 165580141 (95 llamadas)
0 segundos

[amchacon]

¿46 segundos?   A mí me ha tardado 2-3 segundos el método clásico.

Excelente código, aunque voy a ser un poco troll y voy a mejorar el método clásico:

Código (cpp) [Seleccionar] Expandir

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define MIN 37
#define MAX 42

typedef unsigned long long fibo_t;

unsigned long long llamadas_clasico;
unsigned long long llamadas_optimizado;

#define MAX_CACHE 33

unsigned long long Fibonaci_Cache[MAX_CACHE] = {1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,1597
,2584,4181,6765,10946,17711,28657,46368,75025,121393,196418,317811,514229,832040,1346269,2178309,31524578};

fibo_t fibo_clasico(fibo_t n)
{
    llamadas_clasico++;

    if(n < MAX_CACHE)
        return Fibonaci_Cache[n];

    return fibo_clasico(n - 1) + fibo_clasico(n - 2);
}

fibo_t fibo(fibo_t n)
{
   llamadas_optimizado++;

   if(n == 0)
       return 0;
   if(n == 1)
       return 1;
   if(n == 2) //Evitamos recurion infinita si n == 2
       return 1;

   if(n % 2) //n impar y mayor que 2
   {
       fibo_t k1,k2;

       k1 = fibo((n + 1) / 2);
       k2 = fibo((n - 1) / 2);

       return k1 * k1 + k2 * k2;
   }

   //n par y mayor que 2
   return fibo(n / 2) * (fibo((n / 2) + 1) + fibo((n / 2) - 1));
   //si n == 2 fibo(n / 2) * fibo(n / 2 + 1) = fibo(1) * fibo(2)
   //habria recursion infinita si no se pone n == 2 como caso base
}

int main(int argc, char *argv[])
{
   fibo_t i,valor;
   time_t ini;

   printf("Fibonacci clasico:\n");

   ini = time(NULL);

   for(i = MIN ; i < MAX ; i++)
   {
       llamadas_clasico = 0;
       valor = fibo_clasico(i);

       printf("  fibo_clasico(%llu) = %llu (%llu llamadas)\n",i,valor,llamadas_clasico);
   }

   printf("%llu segundos\n",time(NULL) - ini);

   printf("\nFibonacci optimizado:\n");

   ini = time(NULL);

   for(i = MIN ; i < MAX ; i++)
   {
       llamadas_optimizado = 0;
       valor = fibo(i);

       printf("  fibo(%llu) = %llu (%llu llamadas)\n",i,valor,llamadas_optimizado);
   }

   printf("%llu segundos\n",time(NULL) - ini);

   return 0;
}

Ala, ahora los tiempos y las llamadas deberían ser parecidos

[eferion]

Si quieres darle una pequeña vuelta de tuerca más puedes usar los siguientes trucos:

mejor preincrementos a postincrementos:

Código (c) [Seleccionar] Expandir


// esto
++variable;

// es mas optimo que esto ( como norma general, ya que depende de las optimizaciones del compilador )
variable++;[/code ]

Para dividir entre dos, es muuuucho más óptimo usar desplazamientos binarios:

[code=c]
// esto
variable = variable >> 1;

// es más óptimo que esto
variable = variable / 2;

Esto lo podrías optimizar un pelin ( siendo quisquillosos )

Código (c) [Seleccionar] Expandir


// Con código redundante
if(n == 0)
    return 0;
if(n == 1)
    return 1;
if(n == 2) //Evitamos recurion infinita si n == 2
    return 1;

// Optimizado
if(n == 0)
    return 0;
if(n <= 2) // Se evita una comprobacion en cada pasada
    return 1;

Si reutilizas algún cálculo costoso, intenta precalcularlo

Código (c) [Seleccionar] Expandir


// Sin optimizar
return fibo(n / 2) * (fibo((n / 2) + 1) + fibo((n / 2) - 1));

// Optimizado
int temp = n / 2; // o mejor aun temp = n >> 1;
return fibo( temp ) * ( fibo( temp + 1 ) + fibo( temp - 1 ) );

También creo recordar, no estoy 100% seguro, que para saber si un número es divisible entre 2 se puede optimizar también de la siguiente forma ( como no aprovechando desplazamientos binarios ):

Código (c) [Seleccionar] Expandir


// Sin optimizar
if(n % 2) //n impar y mayor que 2

// optimizado
if ( ( n >> 1 << 1 ) == n )

Esto último básicamente consiste en eliminar el bit de menor peso del número ... si es par, ese dígito es cero y el número permanece invariable.

Son truquillos tontos si quieres conseguir arañar hasta la última señal de reloj del procesador.

Eso si, te lo has currado, no voy a ser yo el que te quite el mérito.

Enhorabuena.[/code]

[0xDani]

@eferion, había visto algún artículo sobre optimización de código hasta esos extremos en ensamblador, pero nunca se me había ocurrido hacerlo en C  Supongo que para echar el rato está bien xD

EDIT: Estos tíos sí que lo flipaban con la optimización: http://www.thehackademy.net/madchat/vxdevl/vxmags/ddt-1/DDT%231.2_G

[do-while]

EDIT: Estos tíos sí que lo flipaban con la optimización: http://www.thehackademy.net/madchat/vxdevl/vxmags/ddt-1/DDT%231.2_G

Aunque he leído muy poco, el texto parece entretenido. A ver que se cuentan.

Un amigo está haciendo un trabajo sobre como los compiladores de C (no se si se centra en gcc...) traducen a código máquina y, por lo que me dijo, hay diferencias entre variable++, ++variable,variable = variable + 1 y variable += 1. No se si todas las expresiones anteriores generan distinto código máquina o si alguna es equivalente a alguna de las otras... ¿Sabéis algo mas sobre esto?

¡Saludos!

Acabo de leerlo. No me parece demasiado interesante. Se centra sobre todo en reducir la cantidad de código máquina generado. No busca reducir el coste computacional del algoritmo. Si el código inicial tiene un orden O(n), el final tendrá el mismo, por lo que desde el punto de vista del algoritmo, no gana nada...

[amchacon]

Acabo de leerlo. No me parece demasiado interesante. Se centra sobre todo en reducir la cantidad de código máquina generado. No busca reducir el coste computacional del algoritmo. Si el código inicial tiene un orden O(n), el final tendrá el mismo, por lo que desde el punto de vista del algoritmo, no gana nada...

Todas las optimizaciones en ensamblador son lineales

Un amigo está haciendo un trabajo sobre como los compiladores de C (no se si se centra en gcc...) traducen a código máquina y, por lo que me dijo, hay diferencias entre variable++, ++variable,variable = variable + 1 y variable += 1. No se si todas las expresiones anteriores generan distinto código máquina o si alguna es equivalente a alguna de las otras... ¿Sabéis algo mas sobre esto?

Hay una gran diferencia entre variable = variable +1  y usar los operadores de incremento. La primera expresión tienes que mover la variable en un registro, sumarle 1 y devolverlo a la memoria RAM.

La segunda expresión puede hacerse en una sola sentencia con una sentencia de esamblador (creo que era incr).

De todas formas, si activas las optimización del compilador, probablemente esto te lo haga automático.

[0xDani]

Acabo de leerlo. No me parece demasiado interesante. Se centra sobre todo en reducir la cantidad de código máquina generado. No busca reducir el coste computacional del algoritmo. Si el código inicial tiene un orden O(n), el final tendrá el mismo, por lo que desde el punto de vista del algoritmo, no gana nada...

Claro, yo me refería más bien al tipo de optimizaciones que señalaba @eferion. Y en cuanto a las dudas de los incrementos, en un compilador moderno no tiene importancia, porque van a optimizar el código.

[Shout]

@eferion
Realmente no tiene sentido lo de num >> 1 para dividir entre 2.
Lo digo por esto:

Son 1.000.000.000 iteraciones 

En este benchmark, ">>" ha tardado más, pero seguramente es porque tengo un montón de cosas abiertas, ya que antes no tenía abierto nada (ni siquiera VC++) y ">>" fue un poco más rápido (0.05 segundos), pero aún así creo que no vale la pena complicar el código (legibilidad) por 0.1 segundos. Recuerda, no son 1.000 ni 1.000.000 iteraciones, son ¡1.000.000.000!

[rir3760]

Para dividir entre dos, es muuuucho más óptimo usar desplazamientos binarios:

Código (c) [Seleccionar] Expandir


// esto
variable = variable >> 1;

// es más óptimo que esto
variable = variable / 2;

En el desplazamiento si la variable es de un tipo entero con signo el valor del bit de relleno depende de la implementación, en otras palabras no esta garantizado que el resultado de ambas operaciones (división y desplazamiento) sea el mismo. Si es un tipo sin signo el bit de relleno siempre es cero, ahí no hay problema.

También creo recordar, no estoy 100% seguro, que para saber si un número es divisible entre 2 se puede optimizar también de la siguiente forma ( como no aprovechando desplazamientos binarios ):

Código (c) [Seleccionar] Expandir

// Sin optimizar
if(n % 2) //n impar y mayor que 2

// optimizado
if ( ( n >> 1 << 1 ) == n )

Ya que los números pares (positivos) tienen el bit menos significativo a cero se puede reducir la expresión a "!(n & 1)". Y si invertimos la condición (verificar impares) terminamos utilizando un operador:

Código (c) [Seleccionar] Expandir

if (n & 1){
   /* impares */
}else {
   /* pares */
}

Pero como ya se comento es mejor dejar ciertas mejoras (por ejemplo división vs desplazamiento) al compilador.

Un saludo