src/number.h

   1 #ifdef _WIN32
   2 // min y max entran en conflicto con la windows.h, son rebautizadas en Windows
   3 #define min _cpp_min
   4 #define max _cpp_max
   5 #endif
   6
   7 #include <deque>
   8 #include <utility>
   9 #include <algorithm>
  10 #include <iomanip>
  11 #include <stdint.h>
  12
  13 /* sizeof(E) tiene que ser 2*sizeof(N); y son los tipos nativos con los cuales
  14  * se haran las operaciones mas basicas. */
  15
  16 template < typename N, typename E >
  17 struct number;
  18
  19 template < typename N, typename E >
  20 std::ostream& operator<< (std::ostream& os, const number< N, E >& n);
  21
  22 template < typename N = uint32_t, typename E = uint64_t >
  23 struct number
  24 {
  25
  26         // Tipos
  27         typedef N native_type;
  28         typedef E extended_type;
  29         enum sign_type { positive, negative };
  30         typedef typename std::deque< native_type > chunk_type;
  31         typedef typename chunk_type::size_type size_type;
  32         typedef typename chunk_type::iterator iterator;
  33         typedef typename chunk_type::const_iterator const_iterator;
  34         typedef typename chunk_type::reverse_iterator reverse_iterator;
  35         typedef typename chunk_type::const_reverse_iterator const_reverse_iterator;
  36
  37         // Constructores (después de construído, el chunk siempre tiene al
  38         // menos un elemento).
  39         // Constructor default (1 'átomo con valor 0)
  40         number(): chunk(1, 0) {}
  41
  42         // Constructor a partir de buffer (de 'átomos') y tamaño
  43         // Copia cada elemento del buffer como un 'átomo' del chunk
  44         // (el átomo menos significativo es el chunk[0] == buf[0])
  45         number(native_type* buf, size_type len, sign_type sign = positive):
  46                 chunk(buf, buf + len), sign(sign)
  47         {
  48                 fix_empty();
  49         }
  50
  51         // Constructor a partir de un 'átomo' (lo asigna como único elemento
  52         // del chunk). Copia una vez N en el vector.
  53         number(native_type n, sign_type sign = positive):
  54                 chunk(1, n), sign(sign) {}
  55
  56         // TODO constructor a partir de string.
  57
  58         // Operadores
  59         number& operator++ ()
  60         {
  61                 carry(0);
  62                 return *this;
  63         }
  64
  65         number& operator+= (const number& n);
  66         number& operator*= (const number& n);
  67         number& operator<<= (const size_type n);
  68
  69         // Devuelve referencia a 'átomo' i del chunk (no debería ser necesario
  70         // si la multiplicación es un método de este objeto).
  71         native_type& operator[] (size_type i) {
  72                 return chunk[i];
  73         }
  74
  75         // Iteradores (no deberían ser necesarios)
  76         iterator begin() { return chunk.begin(); }
  77         iterator end() { return chunk.end(); }
  78         const_iterator begin() const { return chunk.begin(); }
  79         const_iterator end() const { return chunk.end(); }
  80         reverse_iterator rbegin() { return chunk.rbegin(); }
  81         reverse_iterator rend() { return chunk.rend(); }
  82         const_reverse_iterator rbegin() const { return chunk.rbegin(); }
  83         const_reverse_iterator rend() const { return chunk.rend(); }
  84
  85         // Friends
  86         template < typename NN, typename EE >
  87         friend std::ostream& operator<< (std::ostream& os, const number< NN, EE>& n);
  88
  89         private:
  90         // Atributos
  91         chunk_type chunk;
  92         sign_type sign;
  93
  94         // Helpers
  95         // Normaliza las longitudes de 2 numbers, completando con 0s a la izquierda
  96         // al más pequeño. Sirve para División y Conquista
  97         number& normalize_length(const number& n);
  98         // parte un número en dos mitades de misma longitud, devuelve un par de
  99         // números con (low, high)
 100         std::pair< number, number > split() const;
 101         // Pone un chunk en 0 para que sea un invariante de representación que
 102         // el chunk no sea vacío (siempre tenga la menos un elemento).
 103         void fix_empty() { if (!chunk.size()) chunk.push_back(0); }
 104         // Propaga carry a partir del 'átomo' i (suma 1 al 'átomo' i propagando
 105         // carry)
 106         void carry(size_type i)
 107         {
 108                 if (chunk.size() > i)
 109                 {
 110                         ++chunk[i];
 111                         if (chunk[i] == 0)
 112                                 carry(i+1); // Overflow
 113                 }
 114                 else
 115                         chunk.push_back(1);
 116         }
 117
 118 };
 119
 120 template < typename N, typename E >
 121 number< N, E >& number< N, E >::operator+= (const number< N, E >& n)
 122 {
 123         native_type c = 0;
 124         size_type ini = 0;
 125         size_type fin = std::min(chunk.size(), n.chunk.size());
 126         size_type i; //problema de VC++, da error de redefinición
 127
 128         // "intersección" entre ambos chunks
 129         // +-----+-----+------+------+
 130         // |     |     |      |      | <--- mio
 131         // +-----+-----+------+------+
 132         // +-----+-----+------+
 133         // |     |     |      |        <--- chunk de n
 134         // +-----+-----+------+
 135         //
 136         // |------------------|
 137         // Esto se procesa en este for
 138         for (i = ini; i < fin; ++i)
 139         {
 140                 chunk[i] += n.chunk[i] + c;
 141                 if ((chunk[i] < n.chunk[i]) || \
 142                                 ( (n.chunk[i] == 0) && c && (chunk[i] == 0) ))
 143                         c = 1; // Overflow
 144                 else
 145                         c = 0; // OK
 146         }
 147
 148         // si mi chunk es más grande que el del otro, sólo me queda
 149         // propagar el carry
 150         if (chunk.size() >= n.chunk.size())
 151         {
 152                 if (c)
 153                         carry(fin); // Propago carry
 154                 return *this;
 155         }
 156
 157         // Hay más
 158         // +-----+-----+------+
 159         // |     |     |      |         <--- mío
 160         // +-----+-----+------+
 161         // +-----+-----+------+------+
 162         // |     |     |      |      |  <--- chunk de n
 163         // +-----+-----+------+------+
 164         //
 165         //                    |------|
 166         //            Esto se procesa en este for
 167         // (suma los chunks de n propagando algún carry si lo había)
 168         ini = fin;
 169         fin = n.chunk.size();
 170         for (i = ini; i < fin; ++i)
 171         {
 172                 chunk.push_back(n.chunk[i] + c); // Agrego nuevo átomo
 173                 if (chunk[i] != 0 || !c)
 174                         c = 0; // OK
 175                 else
 176                         c = 1; // Overflow
 177         }
 178
 179         // Si me queda algún carry colgado, hay que agregar un "átomo"
 180         // más al chunk.
 181         if (c)
 182                 chunk.push_back(1); // Último carry
 183
 184         return *this;
 185 }
 186
 187 template < typename N, typename E >
 188 number< N, E > operator+ (const number< N, E >& n1, const number< N, E >& n2)
 189 {
 190         number< N, E > tmp = n1;
 191         tmp += n2;
 192         return tmp;
 193 }
 194
 195 // efectúa un shifteo a izquierda del chunk, agregando 0s en los casilleros menos significativos
 196 template < typename N, typename E >
 197 number< N, E >& number< N, E >::operator<<= (size_type n)
 198 {
 199         size_type i;
 200         for (i = 0; i < n; i++)
 201         {
 202                 chunk.push_front(0);
 203         }
 204         return *this;
 205 }
 206
 207 template < typename N, typename E >
 208 number< N, E > operator<< (const number< N, E >& n, typename number< N, E >::size_type m)
 209 {
 210         number< N, E > tmp = n;
 211         tmp <<= m;
 212         return tmp;
 213 }
 214
 215 template < typename N, typename E >
 216 std::ostream& operator<< (std::ostream& os, const number< N, E >& n)
 217 {
 218         // FIXME sacar una salida bonita en ASCII =)
 219         for (typename number< N, E >::const_iterator i = n.chunk.begin();
 220                         i != n.chunk.end(); ++i)
 221                 os << std::setfill('0') << std::setw(sizeof(N) * 2) << std::hex
 222                         << *i << " ";
 223         return os;
 224 }
 225
 226 template < typename N, typename E >
 227 number< N, E >& number< N, E >::operator*= (const number< N, E >& n)
 228 {
 229         number < N, E > r_op = n;
 230         normalize_length(n);
 231         n.normalize_length(*this);
 232         *this = divide_n_conquer(*this, n);
 233         return *this;
 234 }
 235
 236 template < typename N, typename E >
 237 number< N, E > operator* (const number< N, E >& n1, const number< N, E >& n2)
 238 {
 239         number< N, E > tmp = n1;
 240         tmp *= n2;
 241         return tmp;
 242 }
 243
 244 template < typename N, typename E >
 245 number< N, E >& number< N, E >::normalize_length(const number< N, E >& n)
 246 {
 247         // si son de distinto tamaño tengo que agregar ceros a la izquierda al
 248         // menor para división y conquista
 249         while (chunk.size() < n.chunk.size())
 250         {
 251                 chunk.push_back(0);
 252         }
 253
 254         // si no tiene cantidad par de números le agrego un atomic_type 0 a la
 255         // izquierda para no tener que contemplar splits de chunks impares
 256         if ((chunk.size() % 2) != 0)
 257         {
 258                 chunk.push_back(0);
 259         }
 260 }
 261
 262 template < typename N, typename E >
 263 std::pair< number< N, E >, number< N, E > > number< N, E >::split() const
 264 {
 265         typedef number< N, E > num_type;
 266         typename num_type::size_type full_size = chunk.size();
 267         typename num_type::size_type halves_size = full_size / 2;
 268         typename num_type::size_type i = 0;
 269
 270         // vacío las mitades
 271         std::pair< num_type, num_type > par;
 272
 273         // la primera mitad va al pedazo inferior
 274         for (i = 0; i < halves_size; i++)
 275         {
 276                 par.first.chunk.push_back(chunk[i]);
 277         }
 278
 279         // la segunda mitad (si full_size es impar es 1 más que la primera
 280         // mitad) va al pedazo superior
 281         for ( ; i < full_size; i++)
 282         {
 283                 par.second.chunk.push_back(chunk[i]);
 284         }
 285         return par;
 286 }
 287
 288 // es el algoritmo de división y conquista, que se llama recursivamente
 289 template < typename N, typename E >
 290 number < N, E > karatsuba(number< N, E > u, number< N, E > v)
 291 {
 292         typedef number< N, E > num_type;
 293
 294         // tomo el chunk size de u (el de v DEBE ser el mismo)
 295         typename num_type::size_type chunk_size = u.chunk.size();
 296
 297         if (chunk_size == 1)
 298         {
 299                 // condición de corte. Ver que por más que tenga 1 único
 300                 // elemento puede "rebalsar" la capacidad del atomic_type,
 301                 // como ser multiplicando 0xff * 0xff usando bytes!!!
 302                 return u.chunk[0] * v.chunk[0];
 303         }
 304
 305         std::pair< num_type, num_type > u12 = u.split();
 306         std::pair< num_type, num_type > v12 = v.split();
 307
 308         // Los nombres M, D y H los puso Rosita en clase, cambiar si se les
 309         // ocurren algunos mejores!
 310         // m = u1*v1
 311         // d = u2*v2
 312         // h = (u1+v1)*(u2+v2) = u1*u2+u1*v2+u2*v1+u2*v2
 313         num_type m = karastuba(u12.first, v12.first);
 314         num_type d = karastuba(u12.second, v12.second);
 315         num_type h = karastuba(u12.first + v12.first,
 316                         u12.second + v12.second);
 317
 318         // H-D-M = u1*u2+u1*v2+u2*v1+u2*v2 - u2*v2 - u1*v1 = u1*v2+u2*v1
 319         // u1*v1 << base^N + u1*v2+u2*v1 << base^N/2 + u2*v2
 320         return (m << chunk_size) + ((h - d - m) << chunk_size / 2) + h;
 321
 322 }
 323