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

## 概要
[乗算編](http://qiita.com/marionette-of-u/items/ec1a47e4b7701de75b63)で用いられたアイディアは[除算にも応用される](http://poset.jp/ompa/node3.html).
被除数 $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` で始端, 終端, 範囲の大きさを保持.

```cpp
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;
};
```

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

```cpp
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;
};
```

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

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

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

```cpp
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` を用いて正規化し, 再帰を終了する.

```cpp
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` へ代入する.

```cpp
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$ を求める.
結果の代入される変数を宣言する.

```cpp
quo_rem qr1;
{
    ...;
}
```

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

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

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

```cpp
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());
```

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

```cpp
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 での全体像.

```cpp
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$ の処理についても同様の流れに基づいて行う.

```cpp
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 で得た結果を結合し, 呼び出し元へ値を返す.

```cpp
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);
```

```cpp
return qr;
```

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

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

実行.

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

呼び出し先は以下.

```cpp
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](https://gist.github.com/marionette-of-u/ad3f26594f48ab7deb61)

## 次回予告
FFT 乗算.