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