/html/espaniol/Apuntes/2014-07-15-MCDBinarioFreePascal_Binary_GCD/MCDBinario-B.jpg

Mcd binario o binary gcd en free pascal.

   Martes 15, de Julio del 2014
 

 

En un artículo anterior explique, como implementar el algoritmo de Euclides para hallar el máximo común divisor, en la explicación anterior se observo que el algoritmo original de Euclides erá muy lento con números bastantes grandes, esto debido a que se hacían muchas restas sucesivas, en vista de ello se busco usar residuos, lo que aceleró bastante bien el algoritmo.

Pero existe un algoritmo conocido como mcd binario, que viene del ingles binary gcd o también conocido como el algoritmo de Stein's. Este algoritmo fue dado a conocer en 1967 por el programador israelí Stein´s, pero al parecer este algoritmo ya era conocido por matemáticos chinos, hace ya bastante tiempo.

Aunque en la teoría este algoritmo resulta ser más rápido que el algoritmo de Euclides implementado con residuos, en la práctica implementarlo resulta ser más difícil, y generalmente se recurre al lenguaje ensamblador para obtener mejores resultados.

Pero veamos en que consiste este algoritmo y el porque resulta ser más lento, si se implementa en un lenguaje de programación de alto nivel como free pascal. Este algoritmo se construye en base a los siguientes teoremas:

  1. Si a,b son pares, entonces el mcd(a,b)=2*mcd(a/2,b/2)
  2. Si a es par y b impar, entonces el mcd(a,b)=mcd(a/2,b)
  3. Si a es impar y b par, entonces el mcd(a,b)=mcd(a,b/2)
  4. Si a y b son impares, entonces el mcd(a,b)=mcd(|a-b|/2,el menor de entre a y b)

El algoritmo consiste en los siguientes pasos, asumiendo que los números a y b son pares entonces usamos el teorema 1, tantas veces hasta que uno de los números sea impar, el número obtenido por la cantidad de veces que se hace esta división, se usará más adelante para multiplicar por una potencia de 2 elevado a ese número obtenido, con el último máximo común divisor obtenido. Si aun a o b es impar usamos el teorema 2 ó el teorema 3, tantas veces hasta que ambos sean impares. Una vez obtenido a y b impares usamos el teorema 4 para luego dependiendo del valor obtenido con |a-b|/2, sea par o impar volvemos a usar uno de los teoremas 2 o 3. El algoritmo termina cuando obtenemos un mcd(1,d)=1 o un mcd(0,d)=d.

Veamos un ejemplo, calcularemos el mcd de 960 y 656, con la ayuda de estos cuatro teoremas.

 

Ejemplo 01

 

Un primer acercamiento a un programa en pseudocódigo podría ser el siguiente:

 

Mientras a y b no sean pares

   a=a/2

   b=b/2

   c=c+1

Fin Mientras

Repetir

   Mientras a no sea par

      a=a/2

   FinMientras

   Mientras b no sea par

      b=b/2

   FinMientras

   Si a>=b

        entonces a=(a-b)/2

        sino b=(b-a)/2

Hasta a<1

mcd=b * (2^c)

 

En este primer pseudocódigo podemos observa que hay muchos bucles mientras, lo que harán que el algoritmo sea más lento. Así que debemos reemplazar esos bucles por algo que nos permita acelerar el algoritmo. Para ello necesitamos adquirir algunos conocimientos previos sobre los números binarios. Para empezar debemos saber como identificar si un número es par o impar, si expresamos un número en binario, podemos observar que cuando el número es par este terminará en cero y cuando es impar terminará en 1. Veamos algunos ejemplos.

 

Ejemplo 01

 

Esto nos puede ayudar a identificar si el número es par o impar, aunque más adelante se verá que no es necesario identificar si el número es par o impar, es necesario saber esto para poder entender lo que viene a continuación.

En el algoritmo hacemos muchas divisiones con 2 cuando el número es par, estas divisiones se pueden hacer más rápido, si desplazamos un bit a la derecha el número que estamos dividiendo.

Veamos un ejemplo, dividamos el número 104 en binario hasta llegar a convertirlo en impar.

 

Ejemplo 01

 

Pero tal como podemos observar hemos echo 3 divisiones con 2, ya que el número en binario 104 tiene 3 ceros por la derecha, es decir de esta observación podemos decir que debemos desplazar bits por la derecha, la cantidad de ceros que tenga este número por la derecha.

Esto lo podemos hacer en freepascal, con el operador shr que nos permite desplazar un bit a la derecha tantas veces queramos. En este ejemplo para convertir 104 a 13, debemos hacer lo siguiente 104 shr 3, usamos el 3 porque 104 tiene 3 ceros por la derecha. Hacer esto equivale a dividir también el número por potencias de 2, es decir 104 shr 3 = 104/(2^3).

Ahora que ya sabemos que podemos usar shr para hacer divisiones con 2 más rápido necesitamos saber cuantos ceros por la derecha tiene un número en binario, para poder dividir el número hasta obtener un número impar. Aquí es en donde necesitamos de una instrucción en ensamblador que nos permita saber cuantos ceros tiene en binario un número por la derecha, en procesadores intel este es conocido como bsf.

En freepascal necesitamos crear la función que nos permita contar estos ceros, la función sería la siguiente:

 

 
function bsf(n:longword):byte;assembler;nostackframe;
  asm
    bsfl n,%eax
    ret		
  end;
 

Código fuente 1: Función bsf en ensamblador.

 

esta función acepta como parámetro el número n, que es de 32bits (un longword es de 4 bytes), esta función tiene la opción nostackframe, esto hace que la función no guarde nada en la pila permitiendo que la función trabaje más rápido. La instrucción en ensamblador:

 

bsfl n,%eax

 

es la encargada de contar los ceros, este contará los ceros que tiene n y los colocará en el registro del microprocesador eax, es en el registro eax en donde se coloca el valor devuelto por la función, la instrucción ret, sirve para indicar que la función debe retornar a la instrucción siguiente desde donde fue invocada. La instrucción bsf tiene al final una letra l, que le indica a la computadora que tendrá que operar con el número más grande que soporte el microprocesador, la letra l viene de long. Ahora podemos cambiar el pseudocódigo anterior por el siguiente:

 

c=abs(bfs(a)-bfs(b))

Repetir

   a=a shr bfs(a)

   b=b shr bfs(b)

   Si a>=b

       entonces a=(a-b) shr 1

       sino b=(b-a) shr 1

Hasta a<1

mcd=b shl c

 

Aquí podemos observar que el primer bucle mientras se ha reemplazado por la siguiente instrucción.

 

c=abs(bfs(a)-bfs(b))

 

Cómo nosotros necesitamos obtener la cantidad de divisiones realizadas con el número 2, hasta que uno de ellos sea impar, o mejor dicho obtener el número con menor cantidad de ceros por la derecha, es que hacemos una resta entre la cantidad de ceros de ambos números pares.

Pero esto se puede abreviar más si usamos una característica del operador or, convirtiendo la instrucción anterior por esta:

 

c=bfs(a or b)

 

al hacer a or b, estaremos combinando los bits de ambos números de tal modo, que este nuevo número obtenido tendrá siempre los bits ceros del menor número par. Veamos un ejemplo, hagamos una operación or con los siguientes números:

 

Ejemplo 01

 

Tal como se puede observar al usar el operador or, combinará los números de tal modo que cuando uno de ellos tenga un bit en 1 y en el otro un 1 o un 0, colocará un 1 en el nuevo número. En caso ambos sean 0 siempre colocará un 0. Este comportamiento siempre obtendrá la cantidad de ceros por la derecha del menor número par.

Otra mejora que podemos hacer al segundo algoritmo es cambiar el bucle Repetir:

 

Repetir

    a=a shr bfs(a)

    b=b shr bfs(b)

    Si a>=b

        entonces a=(a-b) shr 1

        sino b=(b-a) shr 1

Hasta a<1

 

por este otro:

 

a=a shr bfs(a);

Repetir

    b=b shr bfs(b);

    if (a>b)

        entonces

            aux=b;

            b=a;

            a=aux

            b-=a

Hasta b<1;

 

Con este cambio evitamos ejecutar dos veces la función bfs, dentro del bucle repetir y sólo lo hacemos con el número que realmente es par, con la ayuda de un intercambio y haciendo sólo una resta. Con esto tendremos el siguiente seudocodigo definitivo para ya implementarlo en freepascal.

 

c=bfs(a or b)

a=a shr bfs(a);

Repetir

   b=b shr bfs(b);

   if (a>b)

       entonces

           aux=b;

           b=a;

           a=aux

           b-=a

Hasta b<1;

mcd=b shl c

 

El siguiente programa en freepascal muestra la implementación de este algoritmo, y compara este algoritmo con el algoritmo de Euclides con residuos, haciendo dos pruebas. Pero existe un detalle en como debemos compilar el programa para compilarlo debemos usar la opción -O3, que permitirá optimizar el código ejecutable, quitando una serie de código no necesario o redundante que se produce cuando el código fuente es convertido a código ejecutable. Obviamente el algoritmo para el mcd binario si se compila sin la opción -O3 funcionará más lento en ambas pruebas. Entonces aquí el programa con la implementación correspondiente:

 

 
{$codepage utf8}
{$mode objfpc}

//compilar: fpc -O3 mcdBinario.pp

Uses sysutils;

Function mcdEuclidesResiduos(a,b:longword):longword;
var r:longword;
Begin
   if (a=0) then begin result:=b; exit end;
   if (b=0) then begin result:=a; exit end;
   if a<>b then
    Begin
       r:=a;
       while r>0 do
        Begin
          a:=b;
          b:=r;
          r:=a mod b
        end
    End;
   result:=b
End;

function bfs(dw:longword):byte;assembler;nostackframe;
  asm
    bsfl dw,%eax
    ret
  end;

Function mcdBinario(a,b:longword):longword;
var desp:byte;
    aux:longword;
Begin
 if (a=0) then begin result:=b; exit end;
 if (b=0) then begin result:=a; exit end;
 desp:=bfs(a or b);
 a:=a shr bfs(a);
 Repeat
   b:=b shr bfs(b);
   if (a>b) then begin
                  aux:=b;
                  b:=a;
                  a:=aux
                 end;
   b-=a
 until b<1;
 result:=a shl desp

End;

var t1,t2:double;
    i,j,rango:longword;
Begin
 rango:=10000;

 Writeln('Prueba i j  10000');

  t1:=time;
  for i:=0 to rango do
   for J:=0 to rango do
    mcdEuclidesResiduos(i,j);
  t2:=time;
  Writeln('mcdEuclidesResiduos  = ',FormatDateTime('hh:nn:ss:z',t2-t1));

  t1:=time;
  for i:=0 to rango do
   for J:=0 to rango do
    mcdBinario(i,j);
  t2:=time;
  Writeln('mcdBinario           = ',FormatDateTime('hh:nn:ss:z',t2-t1));

  Writeln('Prueba i i+1 10000000');
  rango:=10000000;

  t1:=time;
  for i:=0 to rango do mcdEuclidesResiduos(i,i+1);
  t2:=time;
   Writeln('mcdEuclidesResiduos  = ',FormatDateTime('hh:nn:ss:z',t2-t1));

  t1:=time;
  for i:=0 to rango do mcdBinario(i,i+1);
  t2:=time;
  Writeln('mcdBinario           = ',FormatDateTime('hh:nn:ss:z',t2-t1));

End.
 

Código fuente 2: Programa mcdBinario.pp que compara el mcd binario con el mcd de Euclides.

 

Las dos pruebas consisten en lo siguiente: La primera consiste en obtener el mcd de números que van desde 0 hasta 10000, usando dos bucles para comparar los números, y la segunda prueba consiste en obtener el mcd de números consecutivos de 0 hasta 10000000, en este caso usando sólo un bucle.

 

Última revisión: 29/07/2014.

 

Referencias:

 

http://lemire.me/blog/archives/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/

http://hbfs.wordpress.com/2013/12/10/the-speed-of-gcd/

http://en.wikipedia.org/wiki/Binary_GCD_algorithm

 

Delicious

 

 

 
 

  COMENTARIOS