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 #ifdef DEBUG
   8 #include <iostream>
   9 #endif
  10
  11 #include <deque>
  12 #include <utility>
  13 #include <algorithm>
  14 #include <iomanip>
  15 #include <string>
  16 #include <sstream>
  17 #include <istream>
  18 #include <ostream>
  19 #include <cassert>
  20
  21 #ifdef _WIN32
  22 // VC++ no tiene la stdint.h, se agrega a mano
  23 #include "stdint.h"
  24 #else
  25 #include <stdint.h>
  26 #endif
  27
  28 enum sign_type { positive, negative };
  29
  30
  31 /* sizeof(E) tiene que ser 2*sizeof(N); y son los tipos nativos con los cuales
  32  * se haran las operaciones mas basicas. */
  33
  34 template < typename N, typename E >
  35 struct number;
  36
  37 template < typename N, typename E >
  38 std::ostream& operator<< (std::ostream& os, const number< N, E >& n);
  39
  40 template < typename N = uint32_t, typename E = uint64_t >
  41 struct number
  42 {
  43
  44         // Tipos
  45         typedef N native_type;
  46         typedef E extended_type;
  47         typedef typename std::deque< native_type > chunk_type;
  48         typedef typename chunk_type::size_type size_type;
  49         typedef typename chunk_type::iterator iterator;
  50         typedef typename chunk_type::const_iterator const_iterator;
  51         typedef typename chunk_type::reverse_iterator reverse_iterator;
  52         typedef typename chunk_type::const_reverse_iterator const_reverse_iterator;
  53
  54         // Constructores (después de construído, el chunk siempre tiene al
  55         // menos un elemento).
  56         // Constructor default (1 'átomo con valor 0)
  57         number(): chunk(1, 0), sign(positive) {}
  58
  59         // Constructor a partir de buffer (de 'átomos') y tamaño
  60         // Copia cada elemento del buffer como un 'átomo' del chunk
  61         // (el átomo menos significativo es el chunk[0] == buf[0])
  62         number(native_type* buf, size_type len, sign_type s = positive):
  63                 chunk(buf, buf + len), sign(s)
  64         {
  65                 fix_empty();
  66         }
  67
  68         // Constructor a partir de un 'átomo' (lo asigna como único elemento
  69         // del chunk). Copia una vez N en el vector.
  70         number(native_type n, sign_type s = positive):
  71                 chunk(1, n), sign(s) {}
  72
  73         number(std::string str);
  74
  75         // Operadores
  76         number& operator++ ()
  77         {
  78                 carry(0);
  79                 return *this;
  80         }
  81
  82         number& operator+=  (const number& n);
  83         number& operator*=  (const number& n);
  84         number& operator<<= (const size_type n);
  85         number& operator-=  (const number& n);
  86         bool    operator==  (const number& n) const;
  87         bool    operator<   (const number& n) const;
  88
  89         // Compara si es menor en módulo.
  90         bool menorEnModuloQue(const number& n) const;
  91
  92         // Devuelve referencia a 'átomo' i del chunk (no debería ser necesario
  93         // si la multiplicación es un método de este objeto).
  94         native_type& operator[] (size_type i)
  95         {
  96                 return chunk[i];
  97         }
  98
  99         // Iteradores (no deberían ser necesarios)
 100         iterator begin() { return chunk.begin(); }
 101         iterator end() { return chunk.end(); }
 102         const_iterator begin() const { return chunk.begin(); }
 103         const_iterator end() const { return chunk.end(); }
 104         reverse_iterator rbegin() { return chunk.rbegin(); }
 105         reverse_iterator rend() { return chunk.rend(); }
 106         const_reverse_iterator rbegin() const { return chunk.rbegin(); }
 107         const_reverse_iterator rend() const { return chunk.rend(); }
 108
 109         // Friends
 110         template < typename NN, typename EE >
 111         friend std::ostream& operator<< (std::ostream& os, const number< NN, EE>& n);
 112
 113         // Atributos
 114         //private:
 115         mutable chunk_type chunk;
 116         sign_type sign;
 117
 118         // Helpers
 119         // parte un número en dos mitades de misma longitud, devuelve un par de
 120         // números con (low, high)
 121         std::pair< number, number > split() const;
 122         // Pone un chunk en 0 para que sea un invariante de representación que
 123         // el chunk no sea vacío (siempre tenga la menos un elemento).
 124         void fix_empty() { if (!chunk.size()) chunk.push_back(0); }
 125         // Propaga carry a partir del 'átomo' i (suma 1 al 'átomo' i propagando
 126         // carry)
 127         void carry(size_type i)
 128         {
 129                 if (chunk.size() > i)
 130                 {
 131                         ++chunk[i];
 132                         if (chunk[i] == 0)
 133                                 carry(i+1); // Overflow
 134                 }
 135                 else
 136                         chunk.push_back(1);
 137         }
 138         // Propaga borrow a partir del 'átomo' i (resta 1 al 'átomo' i
 139         // propagando borrow)
 140         void borrow(size_type i)
 141         {
 142                 // para poder pedir prestado debo tener uno a la izquierda
 143                 assert (chunk.size() >= i);
 144
 145                 if (chunk[i] == 0)
 146                 {
 147                         borrow(i+1); // Overflow, pido prestado
 148                         chunk[i] = ~((N)0); //quedo con el valor máximo
 149                 }
 150                 else
 151                 {
 152                         --chunk[i]; //tengo para dar, pero pierdo uno yo
 153                 }
 154         }
 155         // Verifica si es un número par
 156         bool es_impar() const
 157         {
 158                 return chunk[0] & 1; // Bit menos significativo
 159         }
 160         // Divide por 2.
 161         number dividido_dos() const
 162         {
 163                 number n = *this;
 164                 bool lsb = 0; // bit menos significativo
 165                 bool msb = 0; // bit más significativo
 166                 for (typename chunk_type::reverse_iterator i = n.chunk.rbegin();
 167                                 i != n.chunk.rend(); ++i)
 168                 {
 169                         lsb = *i & 1; // bit menos significativo
 170                         *i >>= 1;     // shift
 171                         // seteo bit más significativo de ser necesario
 172                         if (msb)
 173                                 *i |= 1 << (sizeof(native_type) * 8 - 1);
 174                         msb = lsb;
 175                 }
 176                 return n;
 177         }
 178
 179 };
 180
 181 inline unsigned ascii2uint(char c)
 182 {
 183         return c & 0xF;
 184 }
 185
 186 // Convierte pasando el string a forma polinómica y evalúa el polinómio
 187 // utilizando la regla de Horner:
 188 // Polinomio: str[0] * 10^0 ... str[size-1] * 10^(size-1)
 189 // Paso inicial: *this = str[0]; z = 10; i = n-1
 190 // Paso iterativo: *this = *this * z + str[i-1]; i--
 191 template < typename N, typename E >
 192 number< N, E >::number(std::string str):
 193         chunk(1, 0), sign(positive)
 194 {
 195         if (!str.size())
 196                 return; // Si está vacío, no hace nada
 197
 198         number< N, E > diez = 10u;
 199         std::string::size_type i = 0;
 200         std::string::size_type fin = str.size() - 1;
 201         if (str[0] == '-') // Si es negativo, salteo el primer caracter
 202                 ++i;
 203         chunk[0] = ascii2uint(str[i]);
 204         while (i < fin)
 205         {
 206                 *this = *this * diez + number< N, E >(ascii2uint(str[i+1]));
 207                 ++i;
 208         }
 209         if (str[0] == '-') // Si es negativo, le pongo el signo
 210                 sign = negative;
 211 }
 212
 213 template < typename N, typename E >
 214 number< N, E >& number< N, E >::operator+= (const number< N, E >& n)
 215 {
 216         // Si tienen distinto signo, restamos...
 217         if (sign != n.sign)
 218         {
 219                 if (sign == positive) // n es negativo
 220                 {
 221                         number< N, E > tmp = n;
 222                         tmp.sign = positive;
 223                         *this -= tmp;
 224                 }
 225                 else // n es positivo, yo negativo
 226                 {
 227                         sign = positive;
 228                         *this = n - *this;
 229                 }
 230                 return *this;
 231         }
 232
 233         native_type c = 0;
 234         size_type ini = 0;
 235         size_type fin = std::min(chunk.size(), n.chunk.size());
 236         size_type i; //problema de VC++, da error de redefinición
 237
 238         // "intersección" entre ambos chunks
 239         // +-----+-----+------+------+
 240         // |     |     |      |      | <--- mio
 241         // +-----+-----+------+------+
 242         // +-----+-----+------+
 243         // |     |     |      |        <--- chunk de n
 244         // +-----+-----+------+
 245         //
 246         // |------------------|
 247         // Esto se procesa en este for
 248         for (i = ini; i < fin; ++i)
 249         {
 250                 chunk[i] += n.chunk[i] + c;
 251                 if ((chunk[i] < n.chunk[i]) || \
 252                                 ( c && ((n.chunk[i] + c) == 0)) || \
 253                                 ( (n.chunk[i] == 0) && c && (chunk[i] == 0) ))
 254                         c = 1; // Overflow
 255                 else
 256                         c = 0; // OK
 257         }
 258
 259         // si mi chunk es más grande que el del otro, sólo me queda
 260         // propagar el carry
 261         if (chunk.size() >= n.chunk.size())
 262         {
 263                 if (c)
 264                         carry(fin); // Propago carry
 265                 return *this;
 266         }
 267
 268         // Hay más
 269         // +-----+-----+------+
 270         // |     |     |      |         <--- mío
 271         // +-----+-----+------+
 272         // +-----+-----+------+------+
 273         // |     |     |      |      |  <--- chunk de n
 274         // +-----+-----+------+------+
 275         //
 276         //                    |------|
 277         //            Esto se procesa en este for
 278         // (suma los chunks de n propagando algún carry si lo había)
 279         ini = fin;
 280         fin = n.chunk.size();
 281         for (i = ini; i < fin; ++i)
 282         {
 283                 chunk.push_back(n.chunk[i] + c); // Agrego nuevo átomo
 284                 if (chunk[i] != 0 || !c)
 285                         c = 0; // OK
 286                 else
 287                         c = 1; // Overflow
 288         }
 289
 290         // Si me queda algún carry colgado, hay que agregar un "átomo"
 291         // más al chunk.
 292         if (c)
 293                 chunk.push_back(1); // Último carry
 294
 295         return *this;
 296 }
 297
 298 template < typename N, typename E >
 299 number< N, E > operator+ (const number< N, E >& n1, const number< N, E >& n2)
 300 {
 301         number< N, E > tmp = n1;
 302         tmp += n2;
 303         return tmp;
 304 }
 305
 306 template < typename N, typename E >
 307 number< N, E >& number< N, E >::operator-= (const number< N, E >& n)
 308 {
 309         // minuendo - substraendo
 310         number< N, E > minuend;
 311         number< N, E > subtrahend;
 312
 313         // voy a hacer siempre el mayor menos el menor
 314         if (menorEnModuloQue(n))
 315         {
 316                 minuend = n;
 317                 subtrahend = *this;
 318                 //minuendo < sustraendo => resultado negativo
 319                 minuend.sign = negative;
 320         }
 321         else
 322         {
 323                 minuend = *this;
 324                 subtrahend = n;
 325                 //minuendo > sustraendo => resultado positivo
 326                 minuend.sign = positive;
 327         }
 328
 329         size_type ini = 0;
 330         size_type fin = std::min(minuend.chunk.size(), subtrahend.chunk.size());
 331         size_type i; //problema de VC++, da error de redefinición
 332
 333         //estoy seguro de que minuend > subtrahend, con lo cual itero hasta el
 334         //size del menor de los dos. Si el otro es más grande, puede ser que
 335         //esté lleno de 0's pero no puede ser realmente mayor como cifra
 336         for (i = ini; i < fin; ++i)
 337         {
 338                 // si no alcanza para restar pido prestado
 339                 if ((minuend.chunk[i] < subtrahend.chunk[i]))
 340                 {
 341                         // no puedo pedir si soy el más significativo ...
 342                         assert (i != fin);
 343
 344                         // le pido uno al que me sigue
 345                         minuend.borrow(i+1);
 346                 }
 347
 348                 // es como hacer 24-5: el 4 pide prestado al 2 (borrow(i+1)) y
 349                 // después se hace 4 + (9-5) + 1
 350
 351                 minuend.chunk[i] += (~((N)0) - subtrahend.chunk[i]) + 1;
 352         }
 353
 354         //retorno el minuendo ya restado
 355         *this = minuend;
 356         return *this;
 357 }
 358
 359 template < typename N, typename E >
 360 number< N, E > operator- (const number< N, E >& n1, const number< N, E >& n2)
 361 {
 362         number< N, E > tmp = n1;
 363         tmp -= n2;
 364         return tmp;
 365 }
 366
 367
 368 template < typename N, typename E >
 369 bool number< N, E >::operator< (const number< N, E >& n) const
 370 {
 371         if (sign != n.sign)
 372         {
 373                 if (sign == positive) // yo positivo, n negativo
 374                         return false; // yo soy más grande
 375                 else // yo negagivo, n positivo
 376                         return true; // n es más grande
 377         }
 378
 379         if (sign == negative)  // Si comparamos 2 negativos, usamos
 380                 return !menorEnModuloQue(n); // "lógica inversa"
 381         else
 382                 return menorEnModuloQue(n);
 383 }
 384
 385
 386 template < typename N, typename E >
 387 bool number< N, E >::menorEnModuloQue(const number< N, E >& n) const
 388 {
 389         size_type i; //problema de VC++, da error de redefinición
 390
 391         if (chunk.size() > n.chunk.size()) // yo tengo más elementos
 392         {
 393                 // Recorro los bytes más significativos (que tengo sólo yo)
 394                 for (i = n.chunk.size(); i < chunk.size(); ++i)
 395                 {
 396                         if (chunk[i] != 0) // Si tengo algo distinto a 0
 397                         {
 398                                 return false; // Entonces soy más grande
 399                         }
 400                 }
 401         }
 402         else if (chunk.size() < n.chunk.size()) // n tiene más elementos
 403         {
 404                 // Recorro los bytes más significativos (que tiene sólo n)
 405                 for (i = chunk.size(); i < n.chunk.size(); ++i)
 406                 {
 407                         if (chunk[i] != 0) // Si n tiene algo distinto a 0
 408                         {
 409                                 return true; // Entonces soy más chico
 410                         }
 411                 }
 412         }
 413         // sigo con la intersección de ambos
 414         size_type fin = std::min(chunk.size(), n.chunk.size());
 415         i = fin;
 416         while (i != 0) {
 417                 --i;
 418
 419                 if (chunk[i] < n.chunk[i]) // Si es menor
 420                 {
 421                         return true;
 422                 }
 423                 else if (chunk[i] > n.chunk[i]) // Si es mayor
 424                 {
 425                         return false;
 426                 }
 427                 // Si es igual tengo que seguir viendo
 428         }
 429
 430         return false; // Son iguales
 431 }
 432
 433
 434 // efectúa un shifteo a izquierda del chunk, agregando 0s en los casilleros
 435 // menos significativos
 436 template < typename N, typename E >
 437 number< N, E >& number< N, E >::operator<<= (size_type n)
 438 {
 439         size_type i;
 440         for (i = 0; i < n; i++)
 441         {
 442                 chunk.push_front(0);
 443         }
 444         return *this;
 445 }
 446
 447 template < typename N, typename E >
 448 number< N, E > operator<< (const number< N, E >& n, typename number< N, E >::size_type m)
 449 {
 450         number< N, E > tmp = n;
 451         tmp <<= m;
 452         return tmp;
 453 }
 454
 455 // Este es un 'workarround' HORRIBLE, de lo peor que hicimos en nuestras vidas,
 456 // pero realmente no encontramos manera alguna de convertir un número a un
 457 // string decimal que no requiera de divisiones sucesivas. Para no cambiar la
 458 // semántica del programa, decidimos convertir externamente nuestra salida en
 459 // hexadecimal a decimal utilizando un programa externo (en este caso Python
 460 // porque sabemos que está disponible en el laboratorio B).
 461 //
 462 // Estamos realmente avergonzados de haber tenido que llegar a esto, pero no nos
 463 // imaginamos que iba a sernos tan compleja esta conversión. Y nuevamente
 464 // pedimos disculpas.
 465 template < typename NN, typename EE >
 466 std::ostream& operator<< (std::ostream& os, const number< NN, EE >& n)
 467 {
 468         std::string cmd = "python -c 'print ";
 469         if (n.sign == negative)
 470                 cmd += '-';
 471         cmd += "0x" + numberToHex(n) + "'";
 472
 473         char buf[BUFSIZ];
 474         FILE *ptr;
 475
 476         if ((ptr = popen(cmd.c_str(), "r")) != NULL)
 477                 while (fgets(buf, BUFSIZ, ptr) != NULL)
 478                         os << buf;
 479         pclose(ptr);
 480         return os;
 481 }
 482
 483 template < typename N, typename E >
 484 std::istream& operator>> (std::istream& is, number< N, E >& n)
 485 {
 486         std::string str;
 487         is >> str;
 488         unsigned base = 10;
 489         if (is.flags() & std::ios_base::hex) // Si lo piden en hexa
 490                 base = 16;
 491         if (is.flags() & std::ios_base::oct) // Si lo piden en octal
 492                 base = 8;
 493         n = number< N, E >(str, base);
 494         return is;
 495 }
 496
 497 template < typename N, typename E >
 498 std::string numberToHex(const number< N, E >& n)
 499 {
 500         std::ostringstream os;
 501         typename number< N, E >::const_reverse_iterator i = n.chunk.rbegin();
 502         typename number< N, E >::const_reverse_iterator end = n.chunk.rend();
 503         // Salteo ceros
 504         for (; i != end; ++i)
 505                 if (*i != 0)
 506                         break;
 507         if (i != end) // Si no llegué al final, imprimo sin 'leading zeros'
 508         {
 509                 os << std::hex << *i;
 510                 ++i; // y voy al próximo
 511         }
 512         // imprimo el resto con 'leading zeros'
 513         for (; i != end; ++i)
 514                 os << std::setfill('0') << std::setw(sizeof(N) * 2) << std::hex
 515                         << *i;
 516         return os.str();
 517 }
 518
 519 template < typename N, typename E >
 520 std::string numberToHexDebug(const number< N, E >& n)
 521 {
 522         std::ostringstream os;
 523         // FIXME sacar una salida bonita en ASCII =)
 524         if (n.sign == negative)
 525                 os << "-";
 526         typename number< N, E >::const_reverse_iterator i = n.chunk.rbegin();
 527         typename number< N, E >::const_reverse_iterator end = n.chunk.rend();
 528         for (; i != end; ++i)
 529                 os << std::setfill('0') << std::setw(sizeof(N) * 2) << std::hex
 530                         << *i << " ";
 531         return os.str();
 532 }
 533
 534 template < typename N, typename E >
 535 bool number< N, E >::operator==(const number< N, E >& n) const
 536 {
 537         if (sign != n.sign)
 538         {
 539                 return false;
 540         }
 541
 542         size_type ini = 0;
 543         size_type fin = std::min(chunk.size(), n.chunk.size());
 544         size_type i; //problema de VC++, da error de redefinición
 545
 546         // "intersección" entre ambos chunks
 547         // +-----+-----+------+------+
 548         // |     |     |      |      | <--- mio
 549         // +-----+-----+------+------+
 550         // +-----+-----+------+
 551         // |     |     |      |        <--- chunk de n
 552         // +-----+-----+------+
 553         //
 554         // |------------------|
 555         // Esto se procesa en este for
 556         for (i = ini; i < fin; ++i)
 557         {
 558                 if (chunk[i] != n.chunk[i])
 559                 {
 560                         return false;
 561                 }
 562         }
 563
 564         // si mi chunk es más grande que el del otro, sólo me queda
 565         // ver si el resto es cero.
 566         chunk_type const *chunk_grande = 0;
 567         if (chunk.size() > n.chunk.size())
 568         {
 569                 chunk_grande = &chunk;
 570                 fin = chunk.size() - n.chunk.size();
 571         }
 572         else if (chunk.size() < n.chunk.size())
 573         {
 574                 chunk_grande = &n.chunk;
 575                 fin = n.chunk.size() - chunk.size();
 576         }
 577
 578         // Si tienen tamaños distintos, vemos que el resto sea cero.
 579         if (chunk_grande)
 580         {
 581                 for (; i < fin; ++i) // Sigo desde el i que había quedado
 582                 {
 583                         if ((*chunk_grande)[i] != 0)
 584                         {
 585                                 return false;
 586                         }
 587                 }
 588         }
 589         return true; // Son iguales
 590 }
 591
 592 template < typename N, typename E >
 593 number< N, E >& number< N, E >::operator*= (const number< N, E >& n)
 594 {
 595         *this = naif(*this, n);
 596         return *this;
 597 }
 598
 599 template < typename N, typename E >
 600 number< N, E > operator* (const number< N, E >& n1, const number< N, E >& n2)
 601 {
 602         return naif(n1, n2);
 603 }
 604
 605 template < typename N, typename E >
 606 std::pair< number< N, E >, number< N, E > > number< N, E >::split() const
 607 {
 608         assert(chunk.size() > 1);
 609         typedef number< N, E > num_type;
 610         typename num_type::size_type full_size = chunk.size();
 611         typename num_type::size_type halves_size = full_size / 2;
 612         typename num_type::size_type i = 0;
 613
 614         // vacío las mitades
 615         std::pair< num_type, num_type > par;
 616
 617         // la primera mitad va al pedazo inferior
 618         par.first.chunk[0] = chunk[0];
 619         for (i = 1; i < halves_size; i++)
 620         {
 621                 par.first.chunk.push_back(chunk[i]);
 622         }
 623
 624         // la segunda mitad (si full_size es impar es 1 más que la primera
 625         // mitad) va al pedazo superior
 626         par.second.chunk[0] = chunk[i];
 627         for (i++ ; i < full_size; i++)
 628         {
 629                 par.second.chunk.push_back(chunk[i]);
 630         }
 631         return par;
 632 }
 633
 634
 635 // Lleva los tamaños de chunk a la potencia de 2 más cercana, eliminando o
 636 // agregando ceros.
 637 template < typename N, typename E >
 638 void normalize_length(const number< N, E >& u, const number< N, E >& v)
 639 {
 640         typename number< N, E >::size_type max, p, t, pot2, size_u, size_v;
 641
 642         // Busco el primer chunk no nulo de u
 643         for (size_u = u.chunk.size() - 1; size_u != 0; --size_u)
 644                 if (u.chunk[size_u] != 0)
 645                         break;
 646         size_u++;
 647
 648         // Busco el primer chunk no nulo de v
 649         for (size_v = v.chunk.size() - 1; size_v != 0; --size_v)
 650                 if (v.chunk[size_v] != 0)
 651                         break;
 652         size_v++;
 653
 654         max = std::max(size_u, size_v);
 655
 656         /* Buscamos hacer crecer a ambos a la potencia de 2 mas proxima; para
 657          * lo cual la obtenemos y guardamos en p. */
 658         t = max;
 659         p = 0;
 660         while ((1u << p) < max)
 661                 p++;
 662
 663         /* Ahora guardamos en pot2 el tamaño que deben tener. */
 664         pot2 = 1 << p;
 665
 666         /* Y finalmente hacemos crecer los dos numeros agregando 0s hasta
 667          * completar sus tamaños. */
 668         u.chunk.resize(pot2, 0);
 669         v.chunk.resize(pot2, 0);
 670 }
 671
 672
 673 /* Algoritmo "naif" (por no decir "cabeza" o "bruto") de multiplicacion. */
 674 template < typename N, typename E >
 675 number < N, E > naif(const number< N, E > &u, const number< N, E > &v)
 676 {
 677         typedef number< N, E > num_type;
 678
 679         normalize_length(u, v);
 680
 681         /* como acabo de normalizar los tamaños son iguales */
 682         typename num_type::size_type chunk_size = u.chunk.size();
 683
 684         sign_type sign;
 685
 686         if (u.sign == v.sign) {
 687                 sign = positive;
 688         } else {
 689                 sign = negative;
 690         }
 691
 692         if (chunk_size == 1)
 693         {
 694                 /* Si llegamos a multiplicar dos de tamaño 1, lo que hacemos
 695                  * es usar la multiplicacion nativa del tipo N, guardando el
 696                  * resultado en el tipo E (que sabemos es del doble de tamaño
 697                  * de N, ni mas ni menos).
 698                  * Luego, armamos un objeto number usando al resultado como
 699                  * buffer. Si, es feo.
 700                  */
 701                 E tmp;
 702                 tmp = static_cast< E >(u.chunk[0]) * static_cast< E >(v.chunk[0]);
 703                 num_type tnum = num_type(reinterpret_cast< N* >(&tmp), 2, sign);
 704                 return tnum;
 705         }
 706
 707         std::pair< num_type, num_type > u12 = u.split();
 708         std::pair< num_type, num_type > v12 = v.split();
 709
 710         /* m11 = u1*v1
 711          * m12 = u1*v2
 712          * m21 = u2*v1
 713          * m22 = u2*v2
 714          */
 715         num_type m11 = naif(u12.first, v12.first);
 716         num_type m12 = naif(u12.first, v12.second);
 717         num_type m21 = naif(u12.second, v12.first);
 718         num_type m22 = naif(u12.second, v12.second);
 719
 720         /* u*v = (u1*v1) * 2^n + (u1*v2 + u2*v1) * 2^(n/2) + u2*v2
 721          * PERO! Como los numeros estan "al reves" nos queda:
 722          *     = m22 * 2^n + (m12 + m21) * 2^(n/2) + m11
 723          */
 724         num_type res;
 725         res = m22 << chunk_size;
 726         res = res + ((m12 + m21) << (chunk_size / 2));
 727         res = res + m11;
 728         res.sign = sign;
 729         return res;
 730 }
 731
 732
 733 /* Algoritmo de multiplicacion de Karatsuba-Ofman
 734  * Ver los comentarios del algoritmo naif, es practicamente identico salvo en
 735  * los calculos numericos que se especifican debajo.
 736  */
 737 template < typename N, typename E >
 738 number < N, E > karatsuba(const number< N, E > &u, const number< N, E > &v)
 739 {
 740         typedef number< N, E > num_type;
 741
 742         normalize_length(u, v);
 743
 744         typename num_type::size_type chunk_size = u.chunk.size();
 745
 746         sign_type sign;
 747
 748         if (u.sign == v.sign) {
 749                 sign = positive;
 750         } else {
 751                 sign = negative;
 752         }
 753
 754         if (chunk_size == 1) {
 755                 E tmp;
 756                 tmp = static_cast< E >(u.chunk[0]) * static_cast< E >(v.chunk[0]);
 757                 num_type tnum = num_type(reinterpret_cast< N* >(&tmp), 2, sign);
 758                 return tnum;
 759         }
 760
 761         std::pair< num_type, num_type > u12 = u.split();
 762         std::pair< num_type, num_type > v12 = v.split();
 763
 764         /* Aca esta la gracia de toda la cuestion:
 765          * m = u1*v1
 766          * d = u2*v2
 767          * h = (u1+u2)*(v1+v2) = u1*u2+u1*v2+u2*v1+u2*v2
 768          *
 769          * h - d - m = u1*v2+u2*v1
 770          * u1*v1 << base^N  +  u1*v2+u2*v1 << base^(N/2)  +  u2*v2
 771          * m << base^N  +  (h - d - m) << base^(N/2)  +  d
 772         */
 773         num_type m = karatsuba(u12.first, v12.first);
 774         num_type d = karatsuba(u12.second, v12.second);
 775
 776         num_type sumfst = u12.first + u12.second;
 777         num_type sumsnd = v12.first + v12.second;
 778         num_type h = karatsuba(sumfst, sumsnd);
 779
 780         num_type res, tmp;
 781
 782         /* tmp = h - d - m */
 783         normalize_length(h, d);
 784         tmp = h - d;
 785         normalize_length(tmp, m);
 786         tmp = tmp - m;
 787
 788         /* Resultado final */
 789         res = d << chunk_size;
 790         res += tmp << (chunk_size / 2);
 791         res += m;
 792         res.sign = sign;
 793
 794         return res;
 795 }
 796
 797
 798 /* Potenciacion usando multiplicaciones sucesivas.
 799  * Toma dos parametros u y v, devuelve u^v; asume v positivo.
 800  */
 801 template < typename N, typename E >
 802 number < N, E > pot_ko(number< N, E > &u, number< N, E > &v)
 803 {
 804         assert(v.sign == positive);
 805         number< N, E > res, i;
 806
 807         res = u;
 808         res.sign = u.sign;
 809
 810         for (i = 1; i < v; i += 1) {
 811                 res = karatsuba(res, u);
 812         }
 813
 814         return res;
 815 }
 816
 817 /* Potenciacion usando división y conquista.
 818  * Toma dos parametros u y v, devuelve u^v; asume v positivo.
 819  *
 820  * El pseudocódigo del algoritmo es:
 821  * pot(x, y):
 822  *      if y == 1:
 823  *              return x
 824  *      res = pot(x, y/2)
 825  *      res = res * res
 826  *      if y es impar:
 827  *              res = res * x
 828  *      return res
 829  *
 830  * Es O(n) ya que la ecuación es T(n) = T(n/2) + O(1)
 831  *
 832  * El grafo que lo 'representa' (siendo los nodos el exponente y) algo como:
 833  *
 834  *                      1 3
 835  *                   _/  |  \_
 836  *                 _/    |    \_
 837  *                /      |      \
 838  *               6       1       6
 839  *             /   \           /   \
 840  *            /     \         /     \
 841  *           3       3       3       3
 842  *          /|\     /|\     /|\     /|\
 843  *         2 1 2   2 1 2   2 1 2   2 1 2
 844  *        / \ / \ / \ / \ / \ / \ / \ / \
 845  *        1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 846  *
 847  */
 848 template < typename N, typename E >
 849 number< N, E > pot_dyc_n(const number< N, E > &x, const number< N, E > &y)
 850 {
 851         assert(y.sign == positive);
 852         if (y == number< N, E >(1))
 853         {
 854                 return x;
 855         }
 856         number< N, E > res = pot_dyc_n(x, y.dividido_dos());
 857         res = naif(res, res);
 858         if (y.es_impar())
 859         {
 860                 res = naif(res, x); // Multiplico por el x que falta
 861         }
 862         return res;
 863 }
 864
 865 /* Idem que pot_dyc_n(), pero usa karatsuba() para las multiplicaciones. */
 866 template < typename N, typename E >
 867 number< N, E > pot_dyc_k(const number< N, E > &x, const number< N, E > &y)
 868 {
 869         assert(y.sign == positive);
 870         if (y == number< N, E >(1))
 871         {
 872                 return x;
 873         }
 874         number< N, E > res = pot_dyc_k(x, y.dividido_dos());
 875         res = karatsuba(res, res);
 876         if (y.es_impar())
 877         {
 878                 res = karatsuba(res, x);
 879         }
 880         return res;
 881 }
 882