Programa algebraico de Multiplicación Veloz o Algoritmo de Karatsuba

Imaginemos que tenemos 2 números a,b con una cantidad muy grande de dí­gitos. Debido a las limitaciones de nuestro sistema computacional, no es posible multiplicar directamente y obtener el resultado a*b, es decir, un programa simple como

    #include <iostream>
    using namespace std;
    
    int main()
    {
      long double a, b;
      cin >> a;
      cin >> b;
      cout << a*b << endl;
      return 0;
    }
    

no es capaz de hacer la multiplicación, pues recordemos que si las variables a,b son del tipo long double, tenemos que

    -3.36210314311209350626 \cdot 10^{-4932} \leq ab \leq 1.18973149535723176502 \cdot 10^{4932}

pero es claro que lo anterior no necesariamente es cierto, de hecho, basta tomar a y b como cualquiera de las cotas anteriores para que el producto se salga de este rango, lo que provoca que el programa no pueda hacer el producto y no lo intente (generalmente arroja un “error de segmentación”). La forma correcta de hacerlo y de intentar esquivar las limitaciones que tiene un ordenador es tratar los números como arreglos. Por ejemplo, en base 10 podemos hacer una correspondencia entre el número 981 y el arreglo cuya coordenada 2 es el 9, la coordenada 1 es el 8 y la coordenada 0 es el 1 (la correspondencia se puede hacer al reves para que la posición en el arreglo coincida con la potencia de 10 en la descomposición decimal). Por supuesto que lo más óptimo es usar la mayor base posible para así disminuir la dimensión del arreglo que almacena al número, sin embargo, para no complicar mucho más las cosas trabajaremos en este post solo con la base 10. Con lo ya dicho, podemos hacernos una idea de la necesidad de utilizar algoritmos que operen arreglos como si estos fueran grandes números; particularmente, la multiplicación es una operación simple que se complica bastante cuando intentamos esquivar las limitaciones del ordenador, y más aún, cuando intentamos minimizar el tiempo de ejecución del algoritmo. Lo que intentaré explicarles ahora será el funcionamiento de uno de los algoritmos “rápidos” existentes para multiplicar números enormemente grandes. A continuación puntualizamos las ideas que definen al algoritmo de Karatsuba:

  1. Sean a,b dos enteros de n dígitos (n par);
  2. Definimos k = \frac{n}{2};
  3. Divimos a en 2 números: a_0 (que corresponde a los k dígitos de la derecha) y a_1 (que corresponde a los n-k dígitos de la izquierda);
  4. Del mismo modo dividimos b en b_0 y b_1;
  5. Calculamos p_2 = a_1b_1, p_1 = (a_1+a_0)(b_1+b_0) y p_0 = a_0b_0;
  6. Retornamos ab = p_2\cdot 10^{2k}+[p_1-(p_2+p_0)]\cdot 10^k + p_0;

Este es el método básico de descomposición que propone el algoritmo, pero suele ocurrir que los números p_2, p_1, p_0 siguen siendo demasiado grandes como para operarlos, por lo que ellos se deben descomponer de la misma forma, obteniendo ciertos números [p_2]_2, [p_2]_1, [p_2]_0, [p_1]_2, [p_1]_1, [p_1]_0, [p_0]_2, [p_0]_1, [p_0]_0, los cuales podrían tener el mismo problema y sea necesario descomponer nuevamente (la cantidad de valores aumenta exponencialmente)… la idea es dejar el producto ab en función de valores lo suficientemente pequeños como para trabajar con ellos “sin importar mucho la cantidad de operaciones que se hagan” (se supone que aumentar la cantidad de operaciones nos permite llegar resultado final). A continuación les dejo un programa en C++ que implementa el algoritmo de Karatsuba; este programa viene acompañado de otro programa que genera aleatoriamente números gigantezcos (nosotros tardaríamos mucho en escribir 2 números de 1000 cifras manualmente) que me permitió estudiar el rendimiento y el tiempo de ejecución del algoritmo dependiendo de lo grande que son los números a multiplicar (un punto realmente muy interesante, pero que no abordaremos en este post):

    #include <iostream>
    #include <stdlib.h>
    #include <stdio.h>
    #include <time.h>
    using namespace std;
    #define MIN_EXP 1
    
    void karatsuba(unsigned char *a, unsigned char *b, unsigned char *ret, int d);
    void Multiplicacion_usual(unsigned char *a, unsigned char *b, unsigned char *ret, int d);
    void Acarreo(unsigned char *a, int d);
    void Obtener_Num(unsigned char *a, int *largo_a, int MAX_DIGITS);
    void Mostrar(unsigned char *a, int d);
    
    int main(int argc, char **argv)
    {
      int pri, sec, y, D ;
      pri=atoi(argv[1]);
      sec=atoi(argv[2]);
      y = (pri > sec) ? pri : sec;
      for(D=1; D<y; D*=2);
      int MAX_DIGITS=D;
      unsigned char * a = new unsigned char [MAX_DIGITS];
      unsigned char * b = new unsigned char [MAX_DIGITS];
      unsigned char * r = new unsigned char [6*MAX_DIGITS];
      int largo_a, largo_b, d, i ;
      clock_t start;
      clock_t stop;
      Obtener_Num(a, &largo_a, MAX_DIGITS);
      Obtener_Num(b, &largo_b, MAX_DIGITS);
      if(largo_a<0 || largo_b<0)
      {
        printf("0\n");
        exit(0);
        return 0;
      }
      i = (largo_a > largo_b) ? largo_a : largo_b;
      for(d=1; d<i; d*=2);
      for(i=largo_a; i<d; i++)
      {
        a[i]=0;
      }
      for(i=largo_b; i<d; i++)
      {
        b[i]=0;
      }
      start = clock();
      stop = start + CLOCKS_PER_SEC;
      for(i=0; clock()<stop; i++)
      {
        karatsuba(a, b, r, d);
        Acarreo(r, 2*d);
      }
      start = clock() - start;
      Mostrar(r, 2*d);
      printf(" %f [seg] (%d iteraciones)\n", (double) start / CLOCKS_PER_SEC / i, i);
      delete [] a;
      delete [] b;
      delete [] r;
    }
    
    void karatsuba(unsigned char *a, unsigned char *b, unsigned char *ret, int d)
    {
      int i ;
      unsigned char *ar = &a[0];
      unsigned char *al = &a[d/2];
      unsigned char *br = &b[0];
      unsigned char *bl = &b[d/2];
      unsigned char *asum = &ret[d*5];
      unsigned char *bsum = &ret[d*5+d/2];
      unsigned char *x1 = &ret[d*0];
      unsigned char *x2 = &ret[d*1];
      unsigned char *x3 = &ret[d*2];
      if(d <= MIN_EXP)
      {
        Multiplicacion_usual(a, b, ret, d);
        return;
      }
      for(i=0; i<d/2; i++)
      {
        asum[i]=al[i]+ar[i];
        bsum[i]=bl[i]+br[i];
      }
      karatsuba(ar, br, x1, d/2);
      karatsuba(al, bl, x2, d/2);
      karatsuba(asum, bsum, x3, d/2);
      for(i=0; i<d; i++)
      {
        x3[i]=x3[i]-x1[i]-x2[i];
      }
      for(i=0; i<d; i++)
      {
        ret[i+d/2]+=x3[i];
      }
    }
    
    void Multiplicacion_usual(unsigned char *a, unsigned char *b, unsigned char *ret, int d)
    {
      int i, j ;
      for(i=0; i<2*d; i++)
      {
        ret[i]=0;
      }
      for(i=0; i<d; i++)
      {
        for(j=0; j<d; j++)
        {
          ret[i+j]+=a[i]*b[j];
        }
      }
    }
    void Acarreo(unsigned char *a, int d)
    {
      int c=0, i ;
      for(i=0; i<d; i++)
      {
        a[i]+=c;
        if(a[i] < 0)
        {
          c=-(-(a[i]+1)/10+1);
        }
        else
        {
          c = a[i] / 10;
        }
        a[i] -= c * 10;
      }
    }
    
    void Obtener_Num(unsigned char *a, int *largo_a, int MAX_DIGITS)
    {
      int c, i ;
      *largo_a = 0;
      while(true)
      {
        c=getchar();
        if(c=='\n' || c==EOF)
        {
          break;
        }
        if(*largo_a>=MAX_DIGITS)
        {
          while(c!='\n' && c!=EOF)
          {
            c=getchar();
          }
        }
        a[*largo_a]=c-'0';
        ++(*largo_a);
      }
      for(i=0; i*2<*largo_a-1; i++)
      {
        c=a[i];
        a[i]=a[*largo_a-i-1];
        a[*largo_a-i-1]=c;
      }
    }
    
    void Mostrar(unsigned char *a, int d)
    {
      int i;
      for(i=d-1; i>0; i--)
      {
        if(a[i]!=0)
        {
          break;
        }
      }
      for(; i >= 0; i--)
      {
        printf("%d", a[i]);
      }
    }
    

Para quienes no lo sepan, “acarreo” es el nombre que suele darse a la “reserva” (como en la suma) en un algoritmo matemático. Si “karatsuba” es el nombre del ejecutable (luego de compilar y suponiendo que todo se encuentra en la carpeta HOME), la forma de ejecutar el programa sería “./karatsuba N1 N2” (esto en la consola de Linux), donde N1 y N2 son los números que deseamos multiplicar. Si el resultado del producto es demasiado grande, suele guardarse el resultado directamente en un archivo, digamos “producto.dat”, usando la siguiente línea: “./karatsuba N1 N2 > producto.dat “.

Anuncios

Si te gustó el post, por favor dejanos tu comentario!!

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s