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
  12 #ifdef _WIN32
  13 // VC++ no tiene la stdint.h, se agrega a mano
  14 #include "stdint.h"
  15 #else
  16 #include <stdint.h>
  17 #endif
  18
  19 enum sign_type { positive, negative };
  20
  21
  22 /* sizeof(E) tiene que ser 2*sizeof(N); y son los tipos nativos con los cuales
  23  * se haran las operaciones mas basicas. */
  24
  25 template < typename N, typename E >
  26 struct number;
  27
  28 template < typename N, typename E >
  29 std::ostream& operator<< (std::ostream& os, const number< N, E >& n);
  30
  31 template < typename N = uint32_t, typename E = uint64_t >
  32 struct number
  33 {
  34
  35         // Tipos
  36         typedef N native_type;
  37         typedef E extended_type;
  38         typedef typename std::deque< native_type > chunk_type;
  39         typedef typename chunk_type::size_type size_type;
  40         typedef typename chunk_type::iterator iterator;
  41         typedef typename chunk_type::const_iterator const_iterator;
  42         typedef typename chunk_type::reverse_iterator reverse_iterator;
  43         typedef typename chunk_type::const_reverse_iterator const_reverse_iterator;
  44
  45         // Constructores (después de construído, el chunk siempre tiene al
  46         // menos un elemento).
  47         // Constructor default (1 'átomo con valor 0)
  48         number(): chunk(1, 0) {}
  49
  50         // Constructor a partir de buffer (de 'átomos') y tamaño
  51         // Copia cada elemento del buffer como un 'átomo' del chunk
  52         // (el átomo menos significativo es el chunk[0] == buf[0])
  53         number(native_type* buf, size_type len, sign_type sign = positive):
  54                 chunk(buf, buf + len), sign(sign)
  55         {
  56                 fix_empty();
  57         }
  58
  59         // Constructor a partir de un 'átomo' (lo asigna como único elemento
  60         // del chunk). Copia una vez N en el vector.
  61         number(native_type n, sign_type sign = positive):
  62                 chunk(1, n), sign(sign) {}
  63
  64         number(const std::string& str);
  65
  66         // Operadores
  67         number& operator++ ()
  68         {
  69                 carry(0);
  70                 return *this;
  71         }
  72
  73         number& operator+= (const number& n);
  74         number& operator*= (const number& n);
  75         number& operator<<= (const size_type n);
  76         number& operator-= (const number& n);
  77         bool    operator< (const number& n);
  78
  79         // Devuelve referencia a 'átomo' i del chunk (no debería ser necesario
  80         // si la multiplicación es un método de este objeto).
  81         native_type& operator[] (size_type i)
  82         {
  83                 return chunk[i];
  84         }
  85
  86         // Iteradores (no deberían ser necesarios)
  87         iterator begin() { return chunk.begin(); }
  88         iterator end() { return chunk.end(); }
  89         const_iterator begin() const { return chunk.begin(); }
  90         const_iterator end() const { return chunk.end(); }
  91         reverse_iterator rbegin() { return chunk.rbegin(); }
  92         reverse_iterator rend() { return chunk.rend(); }
  93         const_reverse_iterator rbegin() const { return chunk.rbegin(); }
  94         const_reverse_iterator rend() const { return chunk.rend(); }
  95
  96         // Friends
  97         template < typename NN, typename EE >
  98         friend std::ostream& operator<< (std::ostream& os, const number< NN, EE>& n);
  99
 100         // Atributos
 101         //private:
 102         chunk_type chunk;
 103         sign_type sign;
 104
 105         // Helpers
 106         // parte un número en dos mitades de misma longitud, devuelve un par de
 107         // números con (low, high)
 108         std::pair< number, number > split() const;
 109         // Pone un chunk en 0 para que sea un invariante de representación que
 110         // el chunk no sea vacío (siempre tenga la menos un elemento).
 111         void fix_empty() { if (!chunk.size()) chunk.push_back(0); }
 112         // Propaga carry a partir del 'átomo' i (suma 1 al 'átomo' i propagando
 113         // carry)
 114         void carry(size_type i)
 115         {
 116                 if (chunk.size() > i)
 117                 {
 118                         ++chunk[i];
 119                         if (chunk[i] == 0)
 120                                 carry(i+1); // Overflow
 121                 }
 122                 else
 123                         chunk.push_back(1);
 124         }
 125         // Propaga borrow a partir del 'átomo' i (resta 1 al 'átomo' i propagando
 126         // borrow)
 127         void borrow(size_type i)
 128         {
 129                 if (chunk.size() >= i)
 130                 {
 131                         if (chunk[i] == 0)
 132                         {
 133                                 borrow(i+1); // Overflow, pido prestado
 134                                 chunk[i] = ~((N)0); //quedo con el valor máximo
 135                         }
 136                         else
 137                         {
 138                                 --chunk[i]; //tengo para dar, pero pierdo uno yo
 139                         }
 140                 }
 141                 //else ERROR, están haciendo a-b con a>b
 142         }
 143 };
 144
 145 template < typename N, typename E >
 146 number< N, E >::number(const std::string& origen)
 147 {
 148         const N MAX_N = (~( (N)0 ) );
 149         E increment = 0;
 150         E acum = 0;
 151
 152         unsigned length = origen.length();
 153         unsigned number_offset = 0;
 154
 155         while (number_offset<length)
 156         {
 157                 // si encuentro un signo + ó - corto
 158                 if (!isdigit(origen[length-number_offset-1]))
 159                         break;
 160
 161                 increment = (10*number_offset)*(origen[length-number_offset-1]-'0');
 162                 if ((acum + increment) > MAX_N)
 163                 {
 164                         chunk.push_back(acum);
 165                 }
 166
 167         }
 168
 169
 170 }
 171
 172 template < typename N, typename E >
 173 number< N, E >& number< N, E >::operator+= (const number< N, E >& n)
 174 {
 175         native_type c = 0;
 176         size_type ini = 0;
 177         size_type fin = std::min(chunk.size(), n.chunk.size());
 178         size_type i; //problema de VC++, da error de redefinición
 179
 180         // "intersección" entre ambos chunks
 181         // +-----+-----+------+------+
 182         // |     |     |      |      | <--- mio
 183         // +-----+-----+------+------+
 184         // +-----+-----+------+
 185         // |     |     |      |        <--- chunk de n
 186         // +-----+-----+------+
 187         //
 188         // |------------------|
 189         // Esto se procesa en este for
 190         for (i = ini; i < fin; ++i)
 191         {
 192                 chunk[i] += n.chunk[i] + c;
 193                 if ((chunk[i] < n.chunk[i]) || \
 194                                 ( (n.chunk[i] == 0) && c && (chunk[i] == 0) ))
 195                         c = 1; // Overflow
 196                 else
 197                         c = 0; // OK
 198         }
 199
 200         // si mi chunk es más grande que el del otro, sólo me queda
 201         // propagar el carry
 202         if (chunk.size() >= n.chunk.size())
 203         {
 204                 if (c)
 205                         carry(fin); // Propago carry
 206                 return *this;
 207         }
 208
 209         // Hay más
 210         // +-----+-----+------+
 211         // |     |     |      |         <--- mío
 212         // +-----+-----+------+
 213         // +-----+-----+------+------+
 214         // |     |     |      |      |  <--- chunk de n
 215         // +-----+-----+------+------+
 216         //
 217         //                    |------|
 218         //            Esto se procesa en este for
 219         // (suma los chunks de n propagando algún carry si lo había)
 220         ini = fin;
 221         fin = n.chunk.size();
 222         for (i = ini; i < fin; ++i)
 223         {
 224                 chunk.push_back(n.chunk[i] + c); // Agrego nuevo átomo
 225                 if (chunk[i] != 0 || !c)
 226                         c = 0; // OK
 227                 else
 228                         c = 1; // Overflow
 229         }
 230
 231         // Si me queda algún carry colgado, hay que agregar un "átomo"
 232         // más al chunk.
 233         if (c)
 234                 chunk.push_back(1); // Último carry
 235
 236         return *this;
 237 }
 238
 239 template < typename N, typename E >
 240 number< N, E > operator+ (const number< N, E >& n1, const number< N, E >& n2)
 241 {
 242         number< N, E > tmp = n1;
 243         tmp += n2;
 244         return tmp;
 245 }
 246
 247 template < typename N, typename E >
 248 number< N, E >& number< N, E >::operator-= (const number< N, E >& n)
 249 {
 250         // minuendo - substraendo
 251         number< N, E > minuend;
 252         number< N, E > subtrahend;
 253
 254         // voy a hacer siempre el mayor menos el menor
 255         if (*this < n)
 256         {
 257                 minuend = n;
 258                 subtrahend = *this;
 259                 //minuendo < sustraendo => resultado negativo
 260                 minuend.sign = negative;
 261         }
 262         else
 263         {
 264                 minuend = *this;
 265                 subtrahend = n;
 266                 //minuendo > sustraendo => resultado positivo
 267                 minuend.sign = positive;
 268         }
 269
 270         native_type c = 0;
 271         size_type ini = 0;
 272         size_type fin = std::min(minuend.chunk.size(), subtrahend.chunk.size());
 273         size_type i; //problema de VC++, da error de redefinición
 274
 275         //estoy seguro de que minuend > subtrahend, con lo cual itero hasta el size del
 276         //menor de los dos. Si el otro es más grande, puede ser que esté lleno de 0's pero
 277         //no puede ser realmente mayor como cifra
 278         for (i = ini; i < fin; ++i)
 279         {
 280                 // si no alcanza para restar pido prestado
 281                 if ((minuend.chunk[i] < subtrahend.chunk[i]))
 282                 {
 283                         minuend.borrow(i);
 284                 }
 285
 286                 // resto el chunk i-ésimo
 287                 minuend.chunk[i] -= subtrahend.chunk[i];
 288         }
 289
 290         //retorno el minuendo ya restado
 291         *this = minuend;
 292         return *this;
 293 }
 294
 295 template < typename N, typename E >
 296 number< N, E > operator- (const number< N, E >& n1, const number< N, E >& n2)
 297 {
 298         number< N, E > tmp = n1;
 299         tmp -= n2;
 300         return tmp;
 301 }
 302
 303
 304 template < typename N, typename E >
 305 bool number< N, E >::operator< (const number< N, E >& n)
 306 {
 307         number< N, E > n1 = *this;
 308         number< N, E > n2 = n;
 309
 310         // igualo los largos
 311         normalize_length(n1, n2);
 312
 313         // obtengo el largo
 314         size_type length = n1.chunk.size();
 315         size_type i = length - 1;
 316
 317         // me voy fijando desde "la cifra" más significativa si alguno es menor que el otro
 318         // sigo iterando si son iguales hasta recorrer todo el número hasta la parte menos significativa
 319         while (i > 0)
 320         {
 321                 if (n1[i]<n2[i])
 322                         return true;
 323                 if (n1[i]>n2[i])
 324                         return false;
 325
 326                 i--;
 327         }
 328
 329         // si llegué hasta acá es porque son iguales, por lo tanto no es menor estricto
 330         return false;
 331
 332 }
 333
 334 // efectúa un shifteo a izquierda del chunk, agregando 0s en los casilleros menos significativos
 335 template < typename N, typename E >
 336 number< N, E >& number< N, E >::operator<<= (size_type n)
 337 {
 338         size_type i;
 339         for (i = 0; i < n; i++)
 340         {
 341                 chunk.push_front(0);
 342         }
 343         return *this;
 344 }
 345
 346 template < typename N, typename E >
 347 number< N, E > operator<< (const number< N, E >& n, typename number< N, E >::size_type m)
 348 {
 349         number< N, E > tmp = n;
 350         tmp <<= m;
 351         return tmp;
 352 }
 353
 354 template < typename N, typename E >
 355 std::ostream& operator<< (std::ostream& os, const number< N, E >& n)
 356 {
 357         // FIXME sacar una salida bonita en ASCII =)
 358         if (n.sign == positive)
 359                 os << "+ ";
 360         else
 361                 os << "- ";
 362         for (typename number< N, E >::const_reverse_iterator i = n.chunk.rbegin();
 363                         i != n.chunk.rend(); ++i)
 364                 os << std::setfill('0') << std::setw(sizeof(N) * 2) << std::hex
 365                         << *i << " ";
 366         return os;
 367 }
 368
 369 template < typename N, typename E >
 370 number< N, E >& number< N, E >::operator*= (const number< N, E >& n)
 371 {
 372         //number < N, E > r_op = n;
 373         //normalize_length(n);
 374         //n.normalize_length(*this);
 375         *this = naif(*this, n);
 376         return *this;
 377 }
 378
 379 template < typename N, typename E >
 380 number< N, E > operator* (const number< N, E >& n1, const number< N, E >& n2)
 381 {
 382         return naif(n1, n2);
 383 }
 384
 385 template < typename N, typename E >
 386 std::pair< number< N, E >, number< N, E > > number< N, E >::split() const
 387 {
 388         typedef number< N, E > num_type;
 389         typename num_type::size_type full_size = chunk.size();
 390         typename num_type::size_type halves_size = full_size / 2;
 391         typename num_type::size_type i = 0;
 392
 393         // vacío las mitades
 394         std::pair< num_type, num_type > par;
 395
 396         // la primera mitad va al pedazo inferior
 397         par.first.chunk[0] = chunk[0];
 398         for (i = 1; i < halves_size; i++)
 399         {
 400                 par.first.chunk.push_back(chunk[i]);
 401         }
 402
 403         // la segunda mitad (si full_size es impar es 1 más que la primera
 404         // mitad) va al pedazo superior
 405         par.second.chunk[0] = chunk[i];
 406         for (i++ ; i < full_size; i++)
 407         {
 408                 par.second.chunk.push_back(chunk[i]);
 409         }
 410         return par;
 411 }
 412
 413
 414 template < typename N, typename E >
 415 void normalize_length(number< N, E >& u, number< N, E >& v)
 416 {
 417         typedef number< N, E > num_type;
 418         typename num_type::size_type max, p, t, pot2;
 419
 420         max = std::max(u.chunk.size(), v.chunk.size());
 421
 422         /* Buscamos hacer crecer a ambos a la potencia de 2 mas proxima; para
 423          * lo cual la obtenemos y guardamos en p. */
 424         t = max;
 425         p = 0;
 426         while (t != 0) {
 427                 t = t >> 1;
 428                 p++;
 429         }
 430
 431         /* Ahora guardamos en pot2 el tamaño que deben tener. */
 432         pot2 = 1 << p;
 433
 434         /* Y finalmente hacemos crecer los dos numeros agregando 0s hasta
 435          * completar sus tamaños. */
 436         while (u.chunk.size() < pot2)
 437                 u.chunk.push_back(0);
 438
 439         while (v.chunk.size() < pot2)
 440                 v.chunk.push_back(0);
 441
 442         return;
 443 }
 444
 445
 446 /* Algoritmo "naif" (por no decir "cabeza" o "bruto") de multiplicacion. */
 447 template < typename N, typename E >
 448 number < N, E > naif(const number< N, E > &u, const number< N, E > &v)
 449 {
 450         typedef number< N, E > num_type;
 451
 452         // tomo el chunk size de u (el de v DEBE ser el mismo)
 453         typename num_type::size_type chunk_size = u.chunk.size();
 454
 455         sign_type sign;
 456
 457         if ( (u.sign == positive && v.sign == positive) ||
 458                         (u.sign == negative && v.sign == negative) ) {
 459                 sign = positive;
 460         } else {
 461                 sign = negative;
 462         }
 463
 464         //printf("naif %d %d\n", u.chunk.size(), v.chunk.size() );
 465
 466         if (chunk_size == 1)
 467         {
 468                 /* Si llegamos a multiplicar dos de tamaño 1, lo que hacemos
 469                  * es usar la multiplicacion nativa del tipo N, guardando el
 470                  * resultado en el tipo E (que sabemos es del doble de tamaño
 471                  * de N, ni mas ni menos).
 472                  * Luego, armamos un objeto number usando al resultado como
 473                  * buffer. Si, es feo.
 474                  */
 475                 E tmp;
 476                 tmp = static_cast< E >(u.chunk[0]) * static_cast< E >(v.chunk[0]);
 477                 num_type tnum = num_type(reinterpret_cast< N* >(&tmp), 2, sign);
 478                 //std::cout << "T:" << tnum << " " << tmp << "\n";
 479                 //printf("1: %lu %lu %llu\n", u.chunk[0], v.chunk[0], tmp);
 480                 return tnum;
 481         }
 482
 483         std::pair< num_type, num_type > u12 = u.split();
 484         std::pair< num_type, num_type > v12 = v.split();
 485
 486         //std::cout << "u:" << u12.first << " - " << u12.second << "\n";
 487         //std::cout << "v:" << v12.first << " - " << v12.second << "\n";
 488
 489         /* m11 = u1*v1
 490          * m12 = u1*v2
 491          * m21 = u2*v1
 492          * m22 = u2*v2
 493          */
 494         num_type m11 = naif(u12.first, v12.first);
 495         num_type m12 = naif(u12.first, v12.second);
 496         num_type m21 = naif(u12.second, v12.first);
 497         num_type m22 = naif(u12.second, v12.second);
 498
 499         /*
 500         printf("csize: %d\n", chunk_size);
 501         std::cout << "11 " << m11 << "\n";
 502         std::cout << "12 " << m12 << "\n";
 503         std::cout << "21 " << m21 << "\n";
 504         std::cout << "22 " << m22 << "\n";
 505         */
 506
 507         /* u*v = (u1*v1) * 2^n + (u1*v2 + u2*v1) * 2^(n/2) + u2*v2
 508          * PERO! Como los numeros estan "al reves" nos queda:
 509          *     = m22 * 2^n + (m12 + m21) * 2^(n/2) + m11
 510          * FIXME: seria mejor hacer el acomode en la llamada a naif arriba?
 511          */
 512         num_type res;
 513         res = m22 << chunk_size;
 514         res = res + ((m12 + m21) << (chunk_size / 2));
 515         res = res + m11;
 516         res.sign = sign;
 517         /*
 518         std::cout << "r: " << res << "\n";
 519         std::cout << "\n";
 520         */
 521         return res;
 522 }
 523
 524
 525 /* Algoritmo de multiplicacion de Karatsuba-Ofman
 526  * Ver los comentarios del algoritmo naif, es practicamente identico salvo en
 527  * los calculos numericos que se especifican debajo.
 528  */
 529 template < typename N, typename E >
 530 number < N, E > karatsuba(const number< N, E > &u, const number< N, E > &v)
 531 {
 532         typedef number< N, E > num_type;
 533
 534         typename num_type::size_type chunk_size = u.chunk.size();
 535
 536         sign_type sign;
 537
 538         if ( (u.sign == positive && v.sign == positive) ||
 539                         (u.sign == negative && v.sign == negative) ) {
 540                 sign = positive;
 541         } else {
 542                 sign = negative;
 543         }
 544
 545         if (chunk_size == 1) {
 546                 E tmp;
 547                 tmp = static_cast< E >(u.chunk[0]) * static_cast< E >(v.chunk[0]);
 548                 num_type tnum = num_type(static_cast< N* >(&tmp), 2, sign);
 549                 return tnum;
 550         }
 551
 552         std::pair< num_type, num_type > u12 = u.split();
 553         std::pair< num_type, num_type > v12 = v.split();
 554
 555         // Los nombres M, D y H los puso Rosita en clase, cambiar si se les
 556         // ocurren algunos mejores!
 557         // m = u1*v1
 558         // d = u2*v2
 559         // h = (u1+v1)*(u2+v2) = u1*u2+u1*v2+u2*v1+u2*v2
 560         num_type m = karastuba(u12.second, v12.second);
 561         num_type d = karastuba(u12.first, v12.first);
 562         num_type h = karastuba(u12.second + v12.second,
 563                         u12.first + v12.first);
 564
 565         // H-D-M = u1*u2+u1*v2+u2*v1+u2*v2 - u2*v2 - u1*v1 = u1*v2+u2*v1
 566         // u1*v1 << base^N + u1*v2+u2*v1 << base^N/2 + u2*v2
 567         num_type res;
 568         res = (m << chunk_size) + ((h - d - m) << (chunk_size / 2) ) + h;
 569         res.sign = sign;
 570         return res;
 571 }
 572
 573
 574 /* Potenciacion usando multiplicaciones sucesivas.
 575  * Toma dos parametros u y v, devuelve u^v; asume v positivo.
 576  */
 577 template < typename N, typename E >
 578 number < N, E > pot_ko(const number< N, E > &u, const number< N, E > &v)
 579 {
 580         number< N, E > res, i;
 581
 582         res = u;
 583         res.sign = u.sign;
 584
 585         for (i = 1; i < v; i += 1) {
 586                 res *= u;
 587         }
 588
 589         return res;
 590 }
 591