Умножение длинных чисел. Метод Карацубы O(nlog23)
Необходимо умножить два числа большой длины
Идея в том, чтобы заменить одно умножение n-значных чисел тремя умножениями n/2-значных чисел.
(a*10^i + b)(c*10^i + d) = ac*10^2i + (ad + bc)*10^i + bd
(ad + bc) + (ac + bd) = (a + b)(c + d)
(ad + bc) = (a + b)(c + d) - (ac + bd)
Сложность O(nlog23)
Листинг С++
#include <cstdio>
#include <vector>
#include <string>
#include <cmath>
using namespace std;
//--------------------- Defines ------------------------------------------------
#define sint int
#define vsint vector < sint >
#define sz size ()
#define at(v, i) (0 <= (i) && (i) < (v).size () ? v[i] : 0)
#define erase_null(v) while ((v).size () > 1 && (v).back () == 0) (v).pop_back ();
#define next_mod(v, i) if (0 <= (i) + 1 && (i) + 1 < (v).size ()) \
{\
v[(i) + 1] += v[(i)] / 10;\
v[(i)] %= 10; \
}
//--------------------- Shifts, Compare ----------------------------------------
vsint shift_left (vsint s, int k)
{
if (k <= 0) return s;
s.insert (s.begin (), k, 0);
return s;
}
vsint shift_right (vsint s, int k)
{
if (k <= 0) return s;
s.erase (s.begin (), s.begin () + min (k, (int)s.sz));
return s;
}
vsint last (vsint a, int k)
{
return vsint (a.begin (), a.begin () + k - 1);
}
//--------------------- Add, Sub, Mul -------------------------------------
vsint add (vsint &a, vsint &b)
{
vsint c (max (a.sz, b.sz) + 1, 0);
for (int i = 0; i < c.sz - 1; ++ i)
{
c[i] += at (a, i) + at (b, i);
next_mod (c, i);
}
erase_null (c);
return c;
}
vsint sub (vsint &a, vsint &b)
{
vsint c (a.sz, 0);
for (int i = 0; i < c.sz; ++ i)
{
c[i] += at (a, i) - at (b, i);
if (c[i] < 0 && i + 1 < c.sz)
{
c[i] += 10;
-- c[i + 1];
}
}
erase_null (c);
return c;
}
vsint mul (vsint &a, vsint &b)
{
vsint c (a.sz + b.sz + 2, 0);
for (int i = 0; i < a.sz; ++ i)
if (a[i])
for (int j = 0; j < b.sz; ++ j)
if (b[j])
{
c[i + j] += a[i] * b[j];
next_mod (c, i + j);
}
for (int i = 0; i < c.sz; ++ i)
next_mod (c, i);
erase_null (c);
return c;
}
vsint karatsuba_mul (vsint & A, vsint & B)
{
if (A.sz < 100 || B.sz < 100)
return mul (A, B);
int i = (A.sz + B.sz) / 4;
if (i >= A.sz)
i = A.sz / 2;
if (i >= B.sz)
i = B.sz / 2;
vsint a = shift_right (A, i);
vsint b = last (A, i + 1);
vsint c = shift_right (B, i);
vsint d = last (B, i + 1);
vsint ac = karatsuba_mul (a, c);
vsint bd = karatsuba_mul (b, d);
vsint ac_bd = add (ac, bd);
vsint a_b = add (a, b);
vsint c_d = add (c, d);
vsint ab_cd = karatsuba_mul (a_b, c_d);
vsint ad_bc = sub (ab_cd, ac_bd);
ad_bc = shift_left (ad_bc, i);
ac = shift_left (ac, i + i);
ac = add (ac, ad_bc);
return add (ac, bd);
}
17:00
13.07.2010
По всем вопросам обращаться: rumterg@gmail.com