多倍精度整数の基礎 ( Karatsuba除算編 )

概要

乗算編で用いられたアイディアは除算にも応用される.
被除数 $u = u_3 x^3 + u_2 x^2 + u_1 x + u_0$, 除数 $v = v_1 x + v_0$ とある時, 商 $u / v = w = w_1 x + w_0$ において我々にとって関心のある値は $w_1$, $w_0$ で, これら二つを求められれば目的は達成される.

アルゴリズム + コード

Step 0: 宣言など

乗算の時と同様に, 既存の値の表現に対し範囲を用いて計算を行う.
std::vector の const_iterator で始端, 終端, 範囲の大きさを保持.

struct kar_range{
    kar_range() = default;
    kar_range(typename data_type::const_iterator first, typename data_type::const_iterator last) : first(first), last(last), size(last - first){}
    kar_range(typename data_type::const_iterator first, typename data_type::const_iterator last, std::size_t size) : first(first), last(last), size(size){}

    typename data_type::const_iterator first, last;
    std::size_t size;
};

除算結果の商と余りを保持.

struct quo_rem{
    integer quo, rem;

    quo_rem() = default;
    quo_rem(const quo_rem&) = default;
    quo_rem(quo_rem &&other) : quo(std::move(other.quo)), rem(std::move(other.rem)){}
    ~quo_rem() = default;
};

処理を行う関数本体の定義.

static quo_rem kar_rec_div(kar_range lhs, kar_range rhs){ ... }

Step 1: 分割間隔を得る.

除数 $v$ を $v_1$, $v_0$ へそれぞれ分けるために, サイズを 2 の倍数へ繰り上げ, 1/2 にする.

std::size_t n = (rhs.size + (rhs.size % 2)) / 2;

Step 2: 再帰の終了条件

除数が十分に小さい場合, Karatsuba 法ではなく古典的な回復型のアルゴリズムを用いて除算を行う.
関数 unsigned_div は第 1 引数と第 2 引数の範囲を被除数, 第 3 引数と第 4 引数の範囲を除数としこの役割を負う.
但し unsigned_div は処理後に発生する最上位桁から続く 0 を無視するので integer::normalize_data_size を用いて正規化し, 再帰を終了する.

if(n <= 2 || lhs.size > n * 4){
    qr = unsigned_div(lhs.first, lhs.last, rhs.first, rhs.last);
    qr.quo.normalize_data_size();
    qr.rem.normalize_data_size();
    return qr;
}

Step 3: 被除数, 除数を分割する

Step 1 で算出した n に基づいて被除数を n 間隔で 4 つに, 除数を n 間隔で 2 つに分割し, kar_range へ代入する.

kar_range u[4], v[2];
{
    std::size_t a = static_cast<std::size_t>(lhs.last - lhs.first);
    for(std::size_t i = 0; i < 4; ++i){
        std::size_t p = i * n;
        u[i].first = lhs.first + std::min(p, a);
        u[i].last = u[i].first + std::min(n, static_cast<std::size_t>(lhs.last - u[i].first));
        u[i].size = u[i].last - u[i].first;
    }

    a = static_cast<std::size_t>(rhs.last - rhs.first);
    for(std::size_t i = 0; i < 2; ++i){
        std::size_t p = i * n;
        v[i].first = rhs.first + std::min(p, a);
        v[i].last = v[i].first + std::min(n, static_cast<std::size_t>(rhs.last - v[i].first));
        v[i].size = v[i].last - v[i].first;
    }
}

Step 4: 除算を再帰する

$w_1$ を求める.
結果の代入される変数を宣言する.

quo_rem qr1;
{
    ...;
}

$(u_3 x + u_2) / v_1$ を計算する.

kar_range u3_u2(u[2].first, u[3].last);
qr1 = kar_rec_div(u3_u2, v[1]);

得た余り qr1.rem から $v_0 w_1$ を引く.

integer vw;
vw.kar_mul_range(v[0], kar_range(qr1.quo.data.begin(), qr1.quo.data.end()));
qr1.rem.kar_sub(vw.data.begin(), vw.data.end());

引いた結果が負値であれば正値になるまで足す = 除算と乗算を用いて足すべき値を求め, 足す.

if(qr1.rem.sign < 0){
    integer &r1(qr1.rem), prod;
    quo_rem mod = kar_rec_div(kar_range(r1.data.begin(), r1.data.end()), rhs);
    if(!mod.rem.data.empty()){
        mod.quo.unsigned_single_add(1);
    }
    prod.kar_mul_range(kar_range(mod.quo.data.begin(), mod.quo.data.end()), rhs);
    r1.kar_signed_add(prod.data.begin(), prod.data.end());
    qr1.quo.unsigned_sub(mod.quo);
}

この Step での全体像.

quo_rem qr1;
{
    kar_range u3_u2(u[2].first, u[3].last);
    qr1 = kar_rec_div(u3_u2, v[1]);

    qr1.rem.elemental_shift(n);
    qr1.rem.kar_add(u[1].first, u[1].last);

    {
        integer vw;
        vw.kar_mul_range(v[0], kar_range(qr1.quo.data.begin(), qr1.quo.data.end()));
        qr1.rem.kar_sub(vw.data.begin(), vw.data.end());
    }

    if(qr1.rem.sign < 0){
        integer &r1(qr1.rem), prod;
        quo_rem mod = kar_rec_div(kar_range(r1.data.begin(), r1.data.end()), rhs);
        if(!mod.rem.data.empty()){
            mod.quo.unsigned_single_add(1);
        }
        prod.kar_mul_range(kar_range(mod.quo.data.begin(), mod.quo.data.end()), rhs);
        r1.kar_signed_add(prod.data.begin(), prod.data.end());
        qr1.quo.unsigned_sub(mod.quo);
    }
}

$w_0$ の処理についても同様の流れに基づいて行う.

quo_rem qr2;
{
    kar_range r1(qr1.rem.data.begin(), qr1.rem.data.end());
    qr2 = kar_rec_div(r1, v[1]);

    integer &r0(qr2.rem);
    r0.elemental_shift(n);
    r0.kar_add(u[0].first, u[0].last);

    {
        integer vw;
        vw.kar_mul_range(v[0], kar_range(qr2.quo.data.begin(), qr2.quo.data.end()));
        qr2.rem.kar_sub(vw.data.begin(), vw.data.end());
    }

    if(r0.sign < 0){
        integer prod;
        quo_rem mod = kar_rec_div(kar_range(r0.data.begin(), r0.data.end()), rhs);
        if(!mod.rem.data.empty()){
            mod.quo.unsigned_single_add(1);
        }
        prod.kar_mul_range(kar_range(mod.quo.data.begin(), mod.quo.data.end()), rhs);
        r0.kar_signed_add(prod.data.begin(), prod.data.end());
        qr2.quo.unsigned_sub(mod.quo);
    }
}

Step 5: 結合

前回の Step で得た結果を結合し, 呼び出し元へ値を返す.

qr.quo = std::move(qr2.quo);
qr.quo.kar_add(qr1.quo.data.begin(), qr1.quo.data.end(), n);
qr.rem = std::move(qr2.rem);

return qr;

速度計測

以下の環境にて速度の計測を行った.

CPU	MEM	Compiler	OS
Core i7-3770	12GiB	Microsoft Visual C++ 2013	Windows 7

実行.

std::mt19937 random_engine;
std::ofstream ofile("result.txt");
kar_div_time(8192, ofile, random_engine);

呼び出し先は以下.

template<class Out, class RandomEngine>
void kar_div_time(std::size_t num, Out &out, RandomEngine &random_engine){
    integer a, b;
    integer::quo_rem c, d;
    for(std::size_t i = 0; i < num; ++i){
        std::uint16_t n = random_engine() & 0x7FFF;
        a.data.push_back((n << 1) + 1);
        b.data.push_back(n);
    }

    {
        boost::timer t;
        c = integer::div(a, b);
        double u = t.elapsed();
        out << u << std::endl;
    }
    {
        boost::timer t;
        d = integer::modern_div(a, b);
        double u = t.elapsed();
        out << u << std::endl;
    }
    if(c.quo == d.quo && c.rem == d.rem){
        out << "success..." << std::endl;
    }else{
        out << "fail..." << std::endl;
    }
}

結果.
上が回復型, 下が Karatsuba 法による除算となっている.

0.612
0.001
success...

全体像

marionette-of-u / integer.hpp

次回予告

FFT 乗算.