Умножение длинных чисел. Метод Карацубы 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