Renombrar "divide_n_conquer" por "karatsuba".

[z.facultad/75.29/dale.git] / src / number.h

diff --git a/src/number.h b/src/number.h
index b5e433f3b6a53bf60138483991a66603f1b15965..4b77f22c9898b0d5f45635e19f60864fbebb7170 100644 (file)
--- a/src/number.h
+++ b/src/number.h
@@ -1,5 +1,5 @@
 #ifdef _WIN32
 #ifdef _WIN32

-//min y max entran en conflicto con la windows.h, son rebautizadas en Windows
+// min y max entran en conflicto con la windows.h, son rebautizadas en Windows

 #define min _cpp_min
 #define max _cpp_max
 #endif
 #define min _cpp_min
 #define max _cpp_max
 #endif


@@ -8,8 +8,10 @@
 #include <utility>
 #include <algorithm>
 #include <iterator>
 #include <utility>
 #include <algorithm>
 #include <iterator>

+#include <stdint.h>

 
 

-//XXX Pensado para andar con unsigned's (si anda con otra cosa es casualidad =)
+/* sizeof(E) tiene que ser 2*sizeof(N); y son los tipos nativos con los cuales
+ * se haran las operaciones mas basicas. */

 
 template < typename N = uint32_t, typename E = uint64_t >
 struct number
 
 template < typename N = uint32_t, typename E = uint64_t >
 struct number


@@ -26,36 +28,43 @@ struct number
        typedef typename chunk_type::reverse_iterator reverse_iterator;
        typedef typename chunk_type::const_reverse_iterator const_reverse_iterator;
 
        typedef typename chunk_type::reverse_iterator reverse_iterator;
        typedef typename chunk_type::const_reverse_iterator const_reverse_iterator;
 

-       // Constructores (después de construído, el chunk siempre tiene al menos
-       // un elemento).
+       // Constructores (después de construído, el chunk siempre tiene al
+       // menos un elemento).

        // Constructor default (1 'átomo con valor 0)
        number(): chunk(1, 0) {}
        // Constructor default (1 'átomo con valor 0)
        number(): chunk(1, 0) {}

+

        // Constructor a partir de buffer (de 'átomos') y tamaño
        // Copia cada elemento del buffer como un 'átomo' del chunk
        // (el átomo menos significativo es el chunk[0] == buf[0])
        number(native_type* buf, size_type len, sign_type sign = positive):
                chunk(buf, buf + len), sign(sign)
        // Constructor a partir de buffer (de 'átomos') y tamaño
        // Copia cada elemento del buffer como un 'átomo' del chunk
        // (el átomo menos significativo es el chunk[0] == buf[0])
        number(native_type* buf, size_type len, sign_type sign = positive):
                chunk(buf, buf + len), sign(sign)

-               { fix_empty(); }
-       // Constructor a partir de un buffer (de 'átomos') terminado en 0
-       // FIXME (en realidad está 'roto' este contructor porque no puedo
-       // inicializar números con un átomo == 0 en el medio)
-       number(native_type* buf, sign_type sign = positive): sign(sign)
-               { while (*buf) chunk.push_back(*(buf++)); fix_empty(); }
-       // Constructor a partir de un 'átomo' (lo asigna como único elemento del
-       // chunk)
+       {
+               fix_empty();
+       }
+
+       // Constructor a partir de un 'átomo' (lo asigna como único elemento
+       // del chunk). Copia una vez N en el vector.

        number(native_type n, sign_type sign = positive):
        number(native_type n, sign_type sign = positive):

-               chunk(1, n), sign(sign) {} // copia una vez n en el vector
+               chunk(1, n), sign(sign) {}
+

        // TODO constructor a partir de string.
 
        // Operadores
        // TODO constructor a partir de string.
 
        // Operadores

-       number& operator++ () { carry(0); return *this; }
+       number& operator++ ()
+       {
+               carry(0);
+               return *this;
+       }
+

        number& operator+= (const number& n);
        number& operator*= (const number& n);
        number& operator<<= (const size_type n);
 
        // Devuelve referencia a 'átomo' i del chunk (no debería ser necesario
        // si la multiplicación es un método de este objeto).
        number& operator+= (const number& n);
        number& operator*= (const number& n);
        number& operator<<= (const size_type n);
 
        // Devuelve referencia a 'átomo' i del chunk (no debería ser necesario
        // si la multiplicación es un método de este objeto).

-       native_type& operator[] (size_type i) { return chunk[i]; }
+       native_type& operator[] (size_type i) {
+               return chunk[i];
+       }

 
        // Iteradores (no deberían ser necesarios)
        iterator begin() { return chunk.begin(); }
 
        // Iteradores (no deberían ser necesarios)
        iterator begin() { return chunk.begin(); }


@@ -89,9 +98,11 @@ struct number
                if (chunk.size() > i)
                {
                        ++chunk[i];
                if (chunk.size() > i)
                {
                        ++chunk[i];

-                       if (!chunk[i]) carry(i+1); // Overflow
+                       if (chunk[i] == 0)
+                               carry(i+1); // Overflow

                }
                }

-               else chunk.push_back(1);
+               else
+                       chunk.push_back(1);

        }
 
 };
        }
 
 };


@@ -111,22 +122,28 @@ number< N, E >& number< N, E >::operator+= (const number< N, E >& n)
        // +-----+-----+------+
        // |     |     |      |        <--- chunk de n
        // +-----+-----+------+
        // +-----+-----+------+
        // |     |     |      |        <--- chunk de n
        // +-----+-----+------+

-       // 
+       //

        // |------------------|
        // Esto se procesa en este for
        for (i = ini; i < fin; ++i)
        {
                chunk[i] += n.chunk[i] + c;
        // |------------------|
        // Esto se procesa en este for
        for (i = ini; i < fin; ++i)
        {
                chunk[i] += n.chunk[i] + c;

-               if (chunk[i] || (!n.chunk[i] && !c)) c = 0; // OK
-               else                                 c = 1; // Overflow
+               if ((chunk[i] < n.chunk[i]) || \
+                               ( (n.chunk[i] == 0) && c && (chunk[i] == 0) ))
+                       c = 1; // Overflow
+               else
+                       c = 0; // OK

        }
        }

+

        // si mi chunk es más grande que el del otro, sólo me queda
        // propagar el carry
        if (chunk.size() >= n.chunk.size())
        {
        // si mi chunk es más grande que el del otro, sólo me queda
        // propagar el carry
        if (chunk.size() >= n.chunk.size())
        {

-               if (c) carry(fin); // Propago carry
+               if (c)
+                       carry(fin); // Propago carry

                return *this;
        }
                return *this;
        }

+

        // Hay más
        // +-----+-----+------+
        // |     |     |      |         <--- mío
        // Hay más
        // +-----+-----+------+
        // |     |     |      |         <--- mío


@@ -134,7 +151,7 @@ number< N, E >& number< N, E >::operator+= (const number< N, E >& n)
        // +-----+-----+------+------+
        // |     |     |      |      |  <--- chunk de n
        // +-----+-----+------+------+
        // +-----+-----+------+------+
        // |     |     |      |      |  <--- chunk de n
        // +-----+-----+------+------+

-       // 
+       //

        //                    |------|
        //            Esto se procesa en este for
        // (suma los chunks de n propagando algún carry si lo había)
        //                    |------|
        //            Esto se procesa en este for
        // (suma los chunks de n propagando algún carry si lo había)


@@ -143,12 +160,17 @@ number< N, E >& number< N, E >::operator+= (const number< N, E >& n)
        for (i = ini; i < fin; ++i)
        {
                chunk.push_back(n.chunk[i] + c); // Agrego nuevo átomo
        for (i = ini; i < fin; ++i)
        {
                chunk.push_back(n.chunk[i] + c); // Agrego nuevo átomo

-               if (chunk[i] || !c) c = 0; // OK
-               else                c = 1; // Overflow
+               if (chunk[i] != 0 || !c)
+                       c = 0; // OK
+               else
+                       c = 1; // Overflow

        }
        }

+

        // Si me queda algún carry colgado, hay que agregar un "átomo"
        // más al chunk.
        // Si me queda algún carry colgado, hay que agregar un "átomo"
        // más al chunk.

-       if (c) chunk.push_back(1); // Último carry
+       if (c)
+               chunk.push_back(1); // Último carry
+

        return *this;
 }
 
        return *this;
 }
 


@@ -209,14 +231,15 @@ number< N, E > operator* (const number< N, E >& n1, const number< N, E >& n2)
 template < typename N, typename E >
 number< N, E >& number< N, E >::normalize_length(const number< N, E >& n)
 {
 template < typename N, typename E >
 number< N, E >& number< N, E >::normalize_length(const number< N, E >& n)
 {

-       //si son de distinto tamaño tengo que agregar ceros a la izquierda al menor para división y conquista
+       // si son de distinto tamaño tengo que agregar ceros a la izquierda al
+       // menor para división y conquista

        while (chunk.size() < n.chunk.size())
        {
                chunk.push_back(0);
        }
 
        while (chunk.size() < n.chunk.size())
        {
                chunk.push_back(0);
        }
 

-       //si no tiene cantidad par de números le agrego un atomic_type 0 a la izquierda para no tener que contemplar
-       //splits de chunks impares
+       // si no tiene cantidad par de números le agrego un atomic_type 0 a la
+       // izquierda para no tener que contemplar splits de chunks impares

        if ((chunk.size() % 2) != 0)
        {
                chunk.push_back(0);
        if ((chunk.size() % 2) != 0)
        {
                chunk.push_back(0);


@@ -240,7 +263,8 @@ std::pair< number< N, E >, number< N, E > > number< N, E >::split() const
                par.first.chunk.push_back(chunk[i]);
        }
 
                par.first.chunk.push_back(chunk[i]);
        }
 

-       // la segunda mitad (si full_size es impar es 1 más que la primera mitad) va al pedazo superior
+       // la segunda mitad (si full_size es impar es 1 más que la primera
+       // mitad) va al pedazo superior

        for ( ; i < full_size; i++)
        {
                par.second.chunk.push_back(chunk[i]);
        for ( ; i < full_size; i++)
        {
                par.second.chunk.push_back(chunk[i]);


@@ -250,33 +274,37 @@ std::pair< number< N, E >, number< N, E > > number< N, E >::split() const
 
 // es el algoritmo de división y conquista, que se llama recursivamente
 template < typename N, typename E >
 
 // es el algoritmo de división y conquista, que se llama recursivamente
 template < typename N, typename E >

-number < N, E > divide_n_conquer(number< N, E > u, number< N, E > v)
+number < N, E > karatsuba(number< N, E > u, number< N, E > v)

 {
        typedef number< N, E > num_type;
 {
        typedef number< N, E > num_type;

-       //tomo el chunk size de u (el de v DEBE ser el mismo)
+
+       // tomo el chunk size de u (el de v DEBE ser el mismo)

        typename num_type::size_type chunk_size = u.chunk.size();
        typename num_type::size_type chunk_size = u.chunk.size();

-       
+

        if (chunk_size == 1)
        {
        if (chunk_size == 1)
        {

-               //condición de corte. Ver que por más que tenga 1 único elemento puede "rebalsar" la capacidad
-               //del atomic_type, como ser multiplicando 0xff * 0xff usando bytes!!!
+               // condición de corte. Ver que por más que tenga 1 único
+               // elemento puede "rebalsar" la capacidad del atomic_type,
+               // como ser multiplicando 0xff * 0xff usando bytes!!!

                return u.chunk[0] * v.chunk[0];
        }
                return u.chunk[0] * v.chunk[0];
        }

-       
+

        std::pair< num_type, num_type > u12 = u.split();
        std::pair< num_type, num_type > v12 = v.split();
 
        std::pair< num_type, num_type > u12 = u.split();
        std::pair< num_type, num_type > v12 = v.split();
 

-       // Los nombres M, D y H los puso Rosita en clase, cambiar si se les ocurren algunos mejores!
+       // Los nombres M, D y H los puso Rosita en clase, cambiar si se les
+       // ocurren algunos mejores!

        // m = u1*v1
        // d = u2*v2
        // h = (u1+v1)*(u2+v2) = u1*u2+u1*v2+u2*v1+u2*v2
        // m = u1*v1
        // d = u2*v2
        // h = (u1+v1)*(u2+v2) = u1*u2+u1*v2+u2*v1+u2*v2

-       num_type m = divide_n_conquer(u12.first, v12.first);
-       num_type d = divide_n_conquer(u12.second, v12.second);
-       num_type h = divide_n_conquer(u12.first + v12.first, u12.second + v12.second);
+       num_type m = karastuba(u12.first, v12.first);
+       num_type d = karastuba(u12.second, v12.second);
+       num_type h = karastuba(u12.first + v12.first,
+                       u12.second + v12.second);

 
        // H-D-M = u1*u2+u1*v2+u2*v1+u2*v2 - u2*v2 - u1*v1 = u1*v2+u2*v1
        // u1*v1 << base^N + u1*v2+u2*v1 << base^N/2 + u2*v2
        return (m << chunk_size) + ((h - d - m) << chunk_size / 2) + h;
 
        // H-D-M = u1*u2+u1*v2+u2*v1+u2*v2 - u2*v2 - u1*v1 = u1*v2+u2*v1
        // u1*v1 << base^N + u1*v2+u2*v1 << base^N/2 + u2*v2
        return (m << chunk_size) + ((h - d - m) << chunk_size / 2) + h;

-       
+

 }
 
 }