else
chunk.push_back(1);
}
- // Propaga borrow a partir del 'átomo' i (resta 1 al 'átomo' i propagando
- // borrow)
+ // Propaga borrow a partir del 'átomo' i (resta 1 al 'átomo' i
+ // propagando borrow)
void borrow(size_type i)
{
// para poder pedir prestado debo tener uno a la izquierda
assert (chunk.size() >= i);
-
+
if (chunk[i] == 0)
{
borrow(i+1); // Overflow, pido prestado
size_type fin = std::min(minuend.chunk.size(), subtrahend.chunk.size());
size_type i; //problema de VC++, da error de redefinición
- //estoy seguro de que minuend > subtrahend, con lo cual itero hasta el size del
- //menor de los dos. Si el otro es más grande, puede ser que esté lleno de 0's pero
- //no puede ser realmente mayor como cifra
+ //estoy seguro de que minuend > subtrahend, con lo cual itero hasta el
+ //size del menor de los dos. Si el otro es más grande, puede ser que
+ //esté lleno de 0's pero no puede ser realmente mayor como cifra
for (i = ini; i < fin; ++i)
{
// si no alcanza para restar pido prestado
// le pido uno al que me sigue
minuend.borrow(i+1);
}
-
- // es como hacer 24-5: el 4 pide prestado al 2 (borrow(i+1)) y después
- // se hace 4 + (9-5) + 1
- minuend.chunk[i] += (~((N)0) - subtrahend.chunk[i]) + 1;
+ // es como hacer 24-5: el 4 pide prestado al 2 (borrow(i+1)) y
+ // después se hace 4 + (9-5) + 1
+
+ minuend.chunk[i] += (~((N)0) - subtrahend.chunk[i]) + 1;
}
//retorno el minuendo ya restado
return false; // Son iguales
}
-// efectúa un shifteo a izquierda del chunk, agregando 0s en los casilleros menos significativos
+// efectúa un shifteo a izquierda del chunk, agregando 0s en los casilleros
+// menos significativos
template < typename N, typename E >
number< N, E >& number< N, E >::operator<<= (size_type n)
{
chunk_grande = &n.chunk;
fin = n.chunk.size() - chunk.size();
}
- if (chunk_grande) // Si tienen tamaños distintos, vemos que el resto sea cero.
+
+ // Si tienen tamaños distintos, vemos que el resto sea cero.
+ if (chunk_grande)
{
for (; i < fin; ++i) // Sigo desde el i que había quedado
{
}
+// Lleva los tamaños de chunk a la potencia de 2 más cercana, eliminando o
+// agregando ceros.
template < typename N, typename E >
void normalize_length(const number< N, E >& u, const number< N, E >& v)
{
- typedef number< N, E > num_type;
- typename num_type::size_type max, p, t, pot2;
+ typename number< N, E >::size_type max, p, t, pot2, size_u, size_v;
- max = std::max(u.chunk.size(), v.chunk.size());
+ // Busco el primer chunk no nulo de u
+ for (size_u = u.chunk.size() - 1; size_u != 0; --size_u)
+ if (u.chunk[size_u] != 0)
+ break;
+ size_u++;
+
+ // Busco el primer chunk no nulo de v
+ for (size_v = v.chunk.size() - 1; size_v != 0; --size_v)
+ if (v.chunk[size_v] != 0)
+ break;
+ size_v++;
+
+ max = std::max(size_u, size_v);
/* Buscamos hacer crecer a ambos a la potencia de 2 mas proxima; para
* lo cual la obtenemos y guardamos en p. */
t = max;
p = 0;
- while (t != 0) {
- t = t >> 1;
+ while ((1u << p) < max)
p++;
- }
/* Ahora guardamos en pot2 el tamaño que deben tener. */
pot2 = 1 << p;
/* Y finalmente hacemos crecer los dos numeros agregando 0s hasta
* completar sus tamaños. */
- while (u.chunk.size() < pot2)
- u.chunk.push_back(0);
-
- while (v.chunk.size() < pot2)
- v.chunk.push_back(0);
-
- return;
+ u.chunk.resize(pot2, 0);
+ v.chunk.resize(pot2, 0);
}
sign = negative;
}
- //printf("naif %d %d\n", u.chunk.size(), v.chunk.size() );
-
if (chunk_size == 1)
{
/* Si llegamos a multiplicar dos de tamaño 1, lo que hacemos
E tmp;
tmp = static_cast< E >(u.chunk[0]) * static_cast< E >(v.chunk[0]);
num_type tnum = num_type(reinterpret_cast< N* >(&tmp), 2, sign);
- //std::cout << "T:" << tnum << " " << tmp << "\n";
- //printf("1: %lu %lu %llu\n", u.chunk[0], v.chunk[0], tmp);
return tnum;
}
std::pair< num_type, num_type > u12 = u.split();
std::pair< num_type, num_type > v12 = v.split();
- //std::cout << "u:" << u12.first << " - " << u12.second << "\n";
- //std::cout << "v:" << v12.first << " - " << v12.second << "\n";
-
/* m11 = u1*v1
* m12 = u1*v2
* m21 = u2*v1
num_type m21 = naif(u12.second, v12.first);
num_type m22 = naif(u12.second, v12.second);
- /*
- printf("csize: %d\n", chunk_size);
- std::cout << "11 " << m11 << "\n";
- std::cout << "12 " << m12 << "\n";
- std::cout << "21 " << m21 << "\n";
- std::cout << "22 " << m22 << "\n";
- */
-
/* u*v = (u1*v1) * 2^n + (u1*v2 + u2*v1) * 2^(n/2) + u2*v2
* PERO! Como los numeros estan "al reves" nos queda:
* = m22 * 2^n + (m12 + m21) * 2^(n/2) + m11
*/
num_type res;
res = m22 << chunk_size;
- //std::cout << "ra: " << res << "\n";
res = res + ((m12 + m21) << (chunk_size / 2));
- /*
- std::cout << "rb: " << res << "\n";
- std::cout << "12+21: " << (m12 + m21) << "\n";
- std::cout << "cs/2: " << (chunk_size / 2) << "\n";
- std::cout << "t: " << ((m12 + m21) << (chunk_size / 2)) << "\n";
- */
res = res + m11;
- //std::cout << "rc: " << res << "\n";
res.sign = sign;
- //std::cout << "r: " << res << "\n";
- //std::cout << "\n";
return res;
}
std::pair< num_type, num_type > u12 = u.split();
std::pair< num_type, num_type > v12 = v.split();
- /*
- std::cout << "u:" << u12.first << " - " << u12.second << "\n";
- std::cout << "v:" << v12.first << " - " << v12.second << "\n";
- */
-
/* Aca esta la gracia de toda la cuestion:
* m = u1*v1
* d = u2*v2
num_type sumsnd = v12.first + v12.second;
num_type h = karatsuba(sumfst, sumsnd);
- /*
- fflush(stdout); fflush(stderr);
- std::cout << "m: " << m << "\n";
- std::cout << "d: " << d << "\n";
- std::cout << "h: " << h << "\n";
- fflush(stdout); fflush(stderr);
- */
-
num_type res, tmp;
/* tmp = h - d - m */
normalize_length(h, d);
tmp = h - d;
normalize_length(tmp, m);
- /*
- std::cout << "t: " << tmp << "\n";
- std::cout << "m: " << m << "\n";
- */
tmp = tmp - m;
- //std::cout << "t: " << tmp << "\n";
/* Resultado final */
res = d << chunk_size;
number< N, E > pot_dyc_n(const number< N, E > &x, const number< N, E > &y)
{
assert(y.sign == positive);
- //std::cout << "pot(" << x << ", " << y << ")\n";
if (y == number< N, E >(1))
{
- std::cout << "y es 1 => FIN pot(" << x << ", " << y << ")\n";
return x;
}
number< N, E > res = pot_dyc_n(x, y.dividido_dos());
- //std::cout << "y.dividido_dos() = " << y.dividido_dos() << "\n";
- //std::cout << "res = " << res << "\n";
res = naif(res, res);
- //std::cout << "res = " << res << "\n";
if (y.es_impar())
{
- //std::cout << y << " es IMPAR => ";
res = naif(res, x); // Multiplico por el x que falta
- //std::cout << "res = " << res << "\n";
}
- //std::cout << "FIN pot(" << x << ", " << y << ")\n\n";
return res;
}
number< N, E > pot_dyc_k(const number< N, E > &x, const number< N, E > &y)
{
assert(y.sign == positive);
- //std::cout << "pot(" << x << ", " << y << ")\n";
if (y == number< N, E >(1))
{
- //std::cout << "y es 1 => FIN pot(" << x << ", " << y << ")\n";
return x;
}
number< N, E > res = pot_dyc_k(x, y.dividido_dos());
- //std::cout << "y.dividido_dos() = " << y.dividido_dos() << "\n";
- //std::cout << "res = " << res << "\n";
res = karatsuba(res, res);
- //std::cout << "res = " << res << "\n";
if (y.es_impar())
{
- //std::cout << y << " es IMPAR => ";
res = karatsuba(res, x);
- //std::cout << "res = " << res << "\n";
}
- //std::cout << "FIN pot(" << x << ", " << y << ")\n\n";
return res;
}