概要

前回の四則演算編では古典的乗算の実装について触れた. 古典的乗算は基本的に, n 桁の多項式 $f_{n}, g_{n} \in \text{R}[x]$ から得られる積 $h_{2n} = f_{n} g_{n}$ を求めるのに指数時間 = $\text{O}(n^2)$ 回の乗算を必要とする. ここで $x$ は多項式の変数であり $\text{R}$ は乗算と加算において可換で単位元を持つ環 (Ring) とする. 多倍精度整数では $x$ は基数と見なされる.
以下の式は $n$ 次の積 $h_{2n} = f_{n} g_{n}$ を展開した図で 4 回の積が現れ, 更にそれら 4 つの乗算中でも再帰的な分割により, 同様に 4 回の乗算が必要な事が分かる.

f_n g_n = (A_n + B_{n/2})(C_n + D_{n/2}) = A_{2n} C_{2n} + A_n D_n + B_n C_n + B_{n/2} D_{n/2} = h_{2n}

Karatsuba 乗算ではこの $\text{O}(n^2)$ を, 準指数時間 $\text{O}(n^{\text{log}_2 3})$ で計算する.

アルゴリズム

積 $h_{2n} = f_n g_n$ は以下の式に書き変えられる.

f_n g_n = A_{2n} C_{2n} + B_{n/2} D_{n/2} + (A_n + B_n)(C_n + D_n) - A_n C_n - B_n D_n = h_{2n}

加減算は増えるものの, 乗算の回数は $A C$, $B D$, $(A + B)(C + D)$ の 3 回のみになる. そしてこの 3 つの乗算自体も 3 回の乗算に再帰的に処理してしまえば全体の乗算の回数は $\text{O}(n^{\text{log}_2 3})$ になる. これをKaratsuba乗算という.

計測

( 実装の詳細は次の節へ... )
16bit, 128 * (step num) 桁同士の値をそれぞれ古典的乗算 ( 上 ), Karatsuba 乗算 ( 下 ) とで速度比較を行った.
計測を行った環境は処理系 Microsoft Visual C++ 2013 と CPU Core i7-3770 3.40GHz.
速度比較に用いたコードは次の通り.

#include "integer.hpp"

#include <iostream>
#include <boost/timer.hpp>

using integer = multi_precision::integer<>;

void kar_time(std::size_t num){
    integer a, b, c, d;
    for(std::size_t i = 0; i < num; ++i){
        a.data.push_back(i + 1);
        b.data.push_back(num - i);
    }
    {
        boost::timer t;
        c = integer::classic_mul(a, b);
        double u = t.elapsed();
        std::cout << u << std::endl;
    }
    {
        boost::timer t;
        d = integer::mul(a, b);
        double u = t.elapsed();
        std::cout << u << std::endl;
    }
    if(c == d){
        std::cout << "success..." << std::endl;
    }else{
        std::cout << "fail..." << std::endl;
    }
}

int main(){
    for(std::size_t i = 0; i < 16; ++i){
        std::cout << "step " << (i + 1) << "--------" << std::endl;
        kar_time((i + 1) * 128);
        std::cout << std::endl;
    }
    return 0;
}

step 1--------
0.002
0
success...

step 2--------
0.013
0.001
success...

step 3--------
0.038
0.002
success...

step 4--------
0.089
0.002
success...

step 5--------
0.172
0.005
success...

step 6--------
0.295
0.005
success...

step 7--------
0.474
0.005
success...

step 8--------
0.704
0.004
success...

step 9--------
1.002
0.01
success...

step 10--------
1.369
0.01
success...

step 11--------
1.824
0.012
success...

step 12--------
2.363
0.01
success...

step 13--------
2.995
0.013
success...

step 14--------
3.753
0.012
success...

step 15--------
4.608
0.013
success...

step 16--------
5.589
0.008
success...

最終的な 128 * 16 = 2048 桁の乗算では古典的乗算が 5.5 秒以上掛かっているのに対し, Karatsuba乗算は 0.008 秒しか掛かっていない. アルゴリズム内部では $2^n$ 桁の数に対して処理に無駄が無く速いという特徴があるため, 直前の step 15 の 0.013 秒と比べるとしても, 430 倍近くの差が出ている.

準備

range による被演算対象の表現

工夫なしにコードに書き落とした場合, 変数の数が増えデータのコピーなどに時間が掛かり速度が落ちてしまう. これを回避するために変更のない範囲に対して range ベースで各処理を行う方法が有効だと考えられる.

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

2 の冪上への繰り上げ

Karatsuba 乗算では基本的に 2 の冪乗によって range で分割し処理行うため, 単位型に対して繰り上げを行う.
x が 0 になるまで 2 進数での桁数をカウントし, 途中 1 が現れたら n を加算する. n が 1 の時は, 値の最上位以外の bit は 0 のため入力された x を復元して return するのみで良い. n が 2 以上の時は, 値の最上位以外の bit が 1 つ以上 1 であるため繰り上がりをして return する.

static UInt ceil_pow2_single(UInt x){
    std::size_t n = 0, m = 0;
    while(x > 0){
        if((x & 1) > 0){
            ++n;
        }
        ++m;
        x >>= 1;
    }
    if(n == 1){
        return static_cast<UInt>(1 << (m - 1));
    }
    return static_cast<UInt>(1 << m);
}

基数シフト

多倍精度整数値に $\text{Radixs}^n$ を掛ける. これは 2 進数で馴染み深い, 左シフトに相応する. 多倍精度整数では先頭に 0 を n 個追加するのみで良い.

void elemental_shift(std::size_t n){
    data.insert(data.begin(), n, 0);
}

比較関数 kar_less ( range based )

範囲 rhs_first, rhs_last に対して *this < [rhs_first, rhs_last) の比較を行う. kar_sub ( range based ) で使用する.

bool kar_less(typename data_type::const_iterator &rhs_first, typename data_type::const_iterator &rhs_last) const{
    std::size_t lhs_n = data.size(), rhs_n = std::distance(rhs_first, rhs_last);
    typename data_type::const_iterator lhs_first = data.begin(), lhs_iter = lhs_first, rhs_iter = rhs_first;
    if(lhs_n < rhs_n){
        std::size_t n = rhs_n - lhs_n;
        for(std::size_t i = 0; i < n; ++i){
            if(*(rhs_iter + (rhs_n - i - 1)) > 0){
                return true;
            }
        }
        rhs_iter += lhs_n - 1;
        lhs_iter = lhs_first + (lhs_n - 1);
    }else if(lhs_n > rhs_n){
        std::size_t n = lhs_n - rhs_n;
        for(std::size_t i = 0; i < n; ++i){
            if(*(lhs_iter + (lhs_n - i - 1)) > 0){
                return false;
            }
            lhs_iter += rhs_n - 1;
            rhs_iter = rhs_iter + (rhs_n - 1);
        }
    }else{
        lhs_iter += lhs_n - 1;
        rhs_iter += rhs_n - 1;
    }
    for(; ; ){
        UInt l = *lhs_iter, r = *rhs_iter;
        if(l < r){ return true; }
        if(l > r){ return false; }
        if(lhs_iter == lhs_first){
            break;
        }
        --lhs_iter, --rhs_iter;
    }
    return false;
}

減算 kar_sub ( range based )

range ベースの減算処理.

void kar_sub(typename data_type::const_iterator &rhs_first, typename data_type::const_iterator &rhs_last){
    if(kar_less(rhs_first, rhs_last)){
        integer s(*this);
        data.assign(rhs_first, rhs_last);
        unsigned_sub(s);
        sign *=  - 1;
        normalize_data_size();
    }else{
        std::size_t n = std::distance(rhs_first, rhs_last);
        if(n > data.size()){
            unsigned_sub(rhs_first, rhs_first + data.size());
        }else{
            unsigned_sub(rhs_first, rhs_last);
        }
        normalize_data_size();
    }
}

加算 kar_add ( range based )

特に工夫する点はない.

void kar_add(typename data_type::const_iterator &rhs_first, typename data_type::const_iterator &rhs_last){
    std::size_t rhs_size = st::distance(rhs_first, rhs_last);
    if(data.resize() < rhs_size){
        data.resize(rhs_size);
    }
    typename data_type::iterator operand_iter = data.begin();
    DoubleUInt c = 0;
    typename data_type::const_iterator iter = rhs_first, end = rhs_last;
    --end;
    for(; iter != end; ++iter, ++operand_iter){
        UInt &operand(*operand_iter);
        DoubleUInt s = operand + *iter + c;
        operand = static_cast<UInt>(s & base_mask);
        c = s >> BitNum;
    }
    {
        UInt &operand(*operand_iter);
        DoubleUInt s = static_cast<DoubleUInt>(operand) + static_cast<DoubleUInt>(*iter) + c;
        operand = static_cast<UInt>(s & base_mask);
        c = s >> BitNum;
        if(c > 0){
            data.push_back(static_cast<UInt>(c & base_mask));
        }
    }
}

プリミティブな乗算 ( range based )

Karatsuba 乗算で扱えなくなった小さな値同士の乗算. こちらも range ベースなだけの通常の乗算になっている.

void kar_mul_range(const kar_range &lhs, const kar_range &rhs){
    data.resize(lhs.size + rhs.size - 1);
    if(data.empty()){
        normalize_data_size();
        return;
    }
    std::size_t n = 0;
    for(typename data_type::const_iterator rhs_iter = rhs.first; rhs_iter != rhs.last; ++rhs_iter, ++n){
        UInt rhs_value = *rhs_iter, lhs_value;
        std::size_t m = 0;
        for(typename data_type::const_iterator lhs_iter = lhs.first; lhs_iter != lhs.last; ++lhs_iter, ++m){
            lhs_value = *lhs_iter;
            DoubleUInt a = static_cast<DoubleUInt>(rhs_value) * static_cast<DoubleUInt>(lhs_value);
            unsigned_single_add(static_cast<UInt>(a & base_mask), n + m);
            UInt temp = static_cast<UInt>((a >> BitNum) & base_mask);
            if(temp > 0){
                unsigned_single_add(temp, n + m + 1);
            }
        }
    }
}

実装

range の構築と準備の節で書き出した kar_* 系の関数の呼び出しが主な処理となる. 内部で自分自身を 3 回呼び出している点が $\text{O}(n^{log_2 3})$ の根拠になっている.

void kar_rec_mul(const kar_range &x, const kar_range &y){
    std::size_t n = ceil_pow2_single((std::max)(x.size, y.size));
    if(n < 2){
        if(x.size > 0 && y.size > 0){
            kar_mul_range(x, y);
            normalize_data_size();
        }
        return;
    }
    n /= 2;
    std::size_t xn = x.size < n ? x.size : n, yn = y.size < n ? y.size : n;
    kar_range x0, y0, x1, y1;
    x0.first = x.first,       x0.last = x.first + xn,  x0.size = xn;
    x1.first = x.first + xn,  x1.last = x.last,        x1.size = x.size - xn;
    y0.first = y.first,       y0.last = y.first + yn,  y0.size = yn;
    y1.first = y.first + yn,  y1.last = y.last,        y1.size = y.size - yn;
    integer z2, z0;
    {
        integer tx, ty;
        if(x1.first != x1.last){
            tx.data.assign(x1.first, x1.last);
            tx.kar_sub(x0.first, x0.last);
        }else{
            tx.data.assign(x0.first, x0.last);
            tx.sign =  - 1;
        }
        if(y1.first != y1.last){
            ty.data.assign(y1.first, y1.last);
            ty.kar_sub(y0.first, y0.last);
        }else{
            ty.data.assign(y0.first, y0.last);
            ty.sign =  - 1;
        }
        kar_range rx, ry;
        rx.first = tx.data.begin(), rx.last = tx.data.end(), rx.size = tx.data.size();
        ry.first = ty.data.begin(), ry.last = ty.data.end(), ry.size = ty.data.size();
        kar_rec_mul(rx, ry);
        sign = tx.sign != ty.sign ? +1 :  - 1;
    }
    z0.kar_rec_mul(x0, y0);
    add(z0);
    if(x.size >= n && y.size >= n){
        z2.kar_rec_mul(x1, y1);
        add(z2);
    }
    elemental_shift(n);
    add(z0);
    if(x.size >= n && y.size >= n){
        add(z2, n * 2);
    }
    normalize_data_size();
}

全体像

#ifndef MULTI_PRECISION_INCLUDE_INTEGER_HPP
#define MULTI_PRECISION_INCLUDE_INTEGER_HPP

#include <algorithm>
#include <vector>
#include <iostream>
#include <string>
#include <iterator>
#include <utility>
#include <cstdint>
#include <climits>
#include <cmath>

namespace multi_precision{
    template<class UInt = std::uint16_t, class DoubleUInt = std::uint32_t, class DoubleInt = std::int32_t, DoubleUInt BitNum = sizeof(UInt) * CHAR_BIT>
    class integer{
    public:
        static const DoubleUInt base_mask = (static_cast<DoubleUInt>(1) << BitNum) - 1;
        static const UInt half = static_cast<UInt>(static_cast<DoubleUInt>(1) << (BitNum - 1));

        using data_type = std::vector<UInt>;
        data_type data;
        DoubleInt sign = +1;

        integer() = default;
        integer(const integer&) = default;
        integer(integer &&other) : data(std::move(other.data)), sign(other.sign){}
        integer(int x){
            if(x == 0){ return; }
            data.resize(1);
            data[0] = std::abs(x);
            sign = x >= 0 ? +1 :  - 1;
        }

        integer(UInt x){
            if(x == 0){ return; }
            data.resize(1);
            data[0] = x;
            sign = +1;
        }

        integer(const char *str){
            build_from_str(str);
        }

        integer(const std::string &str){
            build_from_str(str.begin());
        }

        template<class StrIter>
        void build_from_str(StrIter iter){
            if(*iter == '-'){
                sign =  - 1;
                ++iter;
            }
            char buff[2] = { 0 };
            while(*iter){
                buff[0] = *iter;
                int a = std::atoi(buff);
                *this *= 10;
                unsigned_single_add(a);
                ++iter;
            }
        }

        ~integer() = default;

        integer &operator =(const integer &other){
            data = other.data;
            sign = other.sign;
            return *this;
        }

        integer &operator =(integer &&other){
            data = std::move(other.data);
            sign = std::move(other.sign);
            return *this;
        }

        bool operator <(const integer &other) const{
            if(sign < other.sign){ return true; }
            if(sign > other.sign){ return false; }
            if(data.size() < other.data.size()){ return sign > 0; }
            if(data.size() > other.data.size()){ return sign < 0; }
            std::size_t i = data.size() - 1;
            do{
                if(data[i] < other.data[i]){ return true; }
            } while(i-- > 0);
            return false;
        }

        bool operator >(const integer &other) const{
            return other < *this;
        }

        bool operator <=(const integer &other) const{
            if(sign < other.sign){ return true; }
            if(sign > other.sign){ return false; }
            if(data.size() < other.data.size()){ return sign > 0; }
            if(data.size() > other.data.size()){ return sign < 0; }
            std::size_t i = data.size() - 1;
            do{
                if(data[i] > other.data[i]){ return false; }
            } while(i-- > 0);
            return true;
        }

        bool operator >=(const integer &other) const{
            return other <= *this;
        }

        bool operator ==(const integer &other) const{
            return sign == other.sign && data == other.data;
        }

        bool operator !=(const integer &other) const{
            return !(*this == other);
        }

        static bool unsigned_less(const data_type &lhs, const data_type &rhs){
            if(lhs.size() < rhs.size()){ return true; }
            if(lhs.size() > rhs.size()){ return false; }
            for(std::size_t i = lhs.size() - 1; i + 1 > 0; --i){
                if(lhs[i] < rhs[i]){ return true; }
                if(lhs[i] > rhs[i]){ return false; }
            }
            return false;
        }

        static bool unsigned_less_eq(const data_type &lhs, const data_type &rhs){
            if(lhs.size() < rhs.size()){ return true; }
            if(lhs.size() > rhs.size()){ return false; }
            for(std::size_t i = lhs.size() - 1; i + 1 > 0; --i){
                if(lhs[i] < rhs[i]){ return true; }
                if(lhs[i] > rhs[i]){ return false; }
            }
            return true;
        }

        void add(const integer &other, std::size_t n = 0){
            if(sign == other.sign){
                unsigned_add(other, n);
            }else{
                if(unsigned_less(data, other.data)){
                    integer temp(*this);
                    data = other.data;
                    unsigned_sub(temp);
                    sign = other.sign;
                }else{
                    unsigned_sub(other);
                }
            }
            normalize_data_size();
        }

        void sub(const integer &other, std::size_t n = 0){
            if(sign != other.sign){
                unsigned_add(other, n);
            }else{
                if(unsigned_less(data, other.data)){
                    integer temp(*this);
                    data = other.data;
                    unsigned_sub(temp);
                    sign = other.sign;
                }else{
                    unsigned_sub(other);
                }
            }
            normalize_data_size();
        }

        void unsigned_add(const integer &other, std::size_t n = 0){
            if(data.size() < other.data.size() + n){
                data.resize(other.data.size() + n);
            }
            DoubleUInt c = 0;
            std::size_t i;
            for(i = 0; i < other.data.size(); ++i){
                DoubleUInt v = static_cast<DoubleUInt>(data[i + n]) + static_cast<DoubleUInt>(other.data[i]) + c;
                data[i + n] = static_cast<UInt>(v & base_mask);
                c = v >> BitNum;
            }
            for(; c > 0; ++i){
                if(i >= data.size()){ data.resize(data.size() + 1); }
                DoubleUInt v = static_cast<DoubleUInt>(data[i + n]) + c;
                data[i + n] = static_cast<UInt>(v & base_mask);
                c = v >> BitNum;
            }
        }

        void unsigned_single_add(UInt c, std::size_t i = 0){
            if(i >= data.size()){ data.resize(i + 1); }
            for(; c > 0; ++i){
                DoubleUInt v = static_cast<DoubleUInt>(data[i]) + c;
                data[i] = static_cast<UInt>(v & base_mask);
                c = v >> BitNum;
            }
        }

        void unsigned_sub(const integer &other, std::size_t i = 0){
            unsigned_sub(other.data.begin(), other.data.end(), i);
        }

        void unsigned_sub(const typename data_type::const_iterator &other_first, const typename data_type::const_iterator &other_last, std::size_t i = 0){
            typename data_type::const_iterator iter = other_first;
            UInt c = 0;
            for(std::size_t n = std::distance(other_first, other_last); i < n; ++i, ++iter){
                const UInt &other_value(*iter);
                UInt t = data[i] - (other_value + c);
                if(data[i] < other_value + c){
                    c = 1;
                }else{
                    c = 0;
                }
                data[i] = t;
            }
            for(; c > 0; ++i){
                UInt t = data[i] - c;
                if(data[i] < c){
                    c = 1;
                }else{
                    c = 0;
                }
                data[i] = t;
            }
            normalize_data_size();
        }

        static void unsigned_mul(integer &r, const integer &lhs, const integer &rhs){
            std::size_t s = lhs.data.size() + rhs.data.size() - 1;
            r.data.resize(s + 1);
            for(std::size_t i = 0; i < lhs.data.size(); ++i){
                for(std::size_t j = 0; j < rhs.data.size(); ++j){
                    DoubleUInt c = static_cast<DoubleUInt>(lhs.data[i]) * static_cast<DoubleUInt>(rhs.data[j]);
                    for(std::size_t k = 0; i + j + k < s + 1; ++k){
                        std::size_t u = i + j + k;
                        DoubleUInt v = static_cast<DoubleUInt>(r.data[u]) + c;
                        UInt a = static_cast<UInt>(v & base_mask);
                        r.data[u] = a;
                        c = v >> BitNum;
                    }
                }
            }
            r.normalize_data_size();
        }
        static integer mul(const integer &lhs, const integer &rhs){
            integer r;
            r.kar_mul(lhs, rhs);
            return r;
        }

        static integer classic_mul(const integer &lhs, const integer &rhs){
            integer r;
            unsigned_mul(r, lhs, rhs);
            return r;
        }

        static void unsigned_single_mul(integer &r, const integer &lhs, UInt rhs){
            std::size_t s = lhs.data.size() + 1;
            r.data.resize(s + 1);
            for(std::size_t i = 0; i < lhs.data.size(); ++i){
                DoubleUInt c = static_cast<DoubleUInt>(lhs.data[i]) * static_cast<DoubleUInt>(rhs);
                for(std::size_t k = 0; i + k < s + 1; ++k){
                    std::size_t u = i + k;
                    DoubleUInt v = static_cast<DoubleUInt>(r.data[u]) + c;
                    UInt a = static_cast<UInt>(v & base_mask);
                    r.data[u] = a;
                    c = v >> BitNum;
                }
            }
            r.normalize_data_size();
        }

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

        void kar_mul(const integer &lhs, const integer &rhs){
            kar_range tx, ty;
            tx.first = lhs.data.begin(), tx.last = lhs.data.end(), tx.size = lhs.data.size();
            ty.first = rhs.data.begin(), ty.last = rhs.data.end(), ty.size = rhs.data.size();
            kar_rec_mul(tx, ty);
            sign = lhs.sign * rhs.sign;
        }


        void kar_mul_range(const kar_range &lhs, const kar_range &rhs){
            data.resize(lhs.size + rhs.size - 1);
            if(data.empty()){
                normalize_data_size();
                return;
            }
            std::size_t n = 0;
            for(typename data_type::const_iterator rhs_iter = rhs.first; rhs_iter != rhs.last; ++rhs_iter, ++n){
                UInt rhs_value = *rhs_iter, lhs_value;
                std::size_t m = 0;
                for(typename data_type::const_iterator lhs_iter = lhs.first; lhs_iter != lhs.last; ++lhs_iter, ++m){
                    lhs_value = *lhs_iter;
                    DoubleUInt a = static_cast<DoubleUInt>(rhs_value) * static_cast<DoubleUInt>(lhs_value);
                    unsigned_single_add(static_cast<UInt>(a & base_mask), n + m);
                    UInt temp = static_cast<UInt>((a >> BitNum) & base_mask);
                    if(temp > 0){
                        unsigned_single_add(temp, n + m + 1);
                    }
                }
            }
        }

        bool kar_less(typename data_type::const_iterator &rhs_first, typename data_type::const_iterator &rhs_last) const{
            std::size_t lhs_n = data.size(), rhs_n = std::distance(rhs_first, rhs_last);
            typename data_type::const_iterator lhs_first = data.begin(), lhs_iter = lhs_first, rhs_iter = rhs_first;
            if(lhs_n < rhs_n){
                std::size_t n = rhs_n - lhs_n;
                for(std::size_t i = 0; i < n; ++i){
                    if(*(rhs_iter + (rhs_n - i - 1)) > 0){
                        return true;
                    }
                }
                rhs_iter += lhs_n - 1;
                lhs_iter = lhs_first + (lhs_n - 1);
            }else if(lhs_n > rhs_n){
                std::size_t n = lhs_n - rhs_n;
                for(std::size_t i = 0; i < n; ++i){
                    if(*(lhs_iter + (lhs_n - i - 1)) > 0){
                        return false;
                    }
                    lhs_iter += rhs_n - 1;
                    rhs_iter = rhs_iter + (rhs_n - 1);
                }
            }else{
                lhs_iter += lhs_n - 1;
                rhs_iter += rhs_n - 1;
            }
            for(; ; ){
                UInt l = *lhs_iter, r = *rhs_iter;
                if(l < r){ return true; }
                if(l > r){ return false; }
                if(lhs_iter == lhs_first){
                    break;
                }
                --lhs_iter, --rhs_iter;
            }
            return false;
        }

        void kar_add(typename data_type::const_iterator &rhs_first, typename data_type::const_iterator &rhs_last){
            std::size_t rhs_size = st::distance(rhs_first, rhs_last);
            if(data.resize() < rhs_size){
                data.resize(rhs_size);
            }
            typename data_type::iterator operand_iter = data.begin();
            DoubleUInt c = 0;
            typename data_type::const_iterator iter = rhs_first, end = rhs_last;
            --end;
            for(; iter != end; ++iter, ++operand_iter){
                UInt &operand(*operand_iter);
                DoubleUInt s = operand + *iter + c;
                operand = static_cast<UInt>(s & base_mask);
                c = s >> BitNum;
            }
            {
                UInt &operand(*operand_iter);
                DoubleUInt s = static_cast<DoubleUInt>(operand) + static_cast<DoubleUInt>(*iter) + c;
                operand = static_cast<UInt>(s & base_mask);
                c = s >> BitNum;
                if(c > 0){
                    data.push_back(static_cast<UInt>(c & base_mask));
                }
            }
        }

        void kar_sub(typename data_type::const_iterator &rhs_first, typename data_type::const_iterator &rhs_last){
            if(kar_less(rhs_first, rhs_last)){
                integer s(*this);
                data.assign(rhs_first, rhs_last);
                unsigned_sub(s);
                sign *=  - 1;
                normalize_data_size();
            }else{
                std::size_t n = std::distance(rhs_first, rhs_last);
                if(n > data.size()){
                    unsigned_sub(rhs_first, rhs_first + data.size());
                }else{
                    unsigned_sub(rhs_first, rhs_last);
                }
                normalize_data_size();
            }
        }

        void kar_rec_mul(const kar_range &x, const kar_range &y){
            std::size_t n = ceil_pow2_single((std::max)(x.size, y.size));
            if(n < 2){
                if(x.size > 0 && y.size > 0){
                    kar_mul_range(x, y);
                    normalize_data_size();
                }
                return;
            }
            n /= 2;
            std::size_t xn = x.size < n ? x.size : n, yn = y.size < n ? y.size : n;
            kar_range x0, y0, x1, y1;
            x0.first = x.first,       x0.last = x.first + xn,  x0.size = xn;
            x1.first = x.first + xn,  x1.last = x.last,        x1.size = x.size - xn;
            y0.first = y.first,       y0.last = y.first + yn,  y0.size = yn;
            y1.first = y.first + yn,  y1.last = y.last,        y1.size = y.size - yn;
            integer z2, z0;
            {
                integer tx, ty;
                if(x1.first != x1.last){
                    tx.data.assign(x1.first, x1.last);
                    tx.kar_sub(x0.first, x0.last);
                }else{
                    tx.data.assign(x0.first, x0.last);
                    tx.sign =  - 1;
                }
                if(y1.first != y1.last){
                    ty.data.assign(y1.first, y1.last);
                    ty.kar_sub(y0.first, y0.last);
                }else{
                    ty.data.assign(y0.first, y0.last);
                    ty.sign =  - 1;
                }
                kar_range rx, ry;
                rx.first = tx.data.begin(), rx.last = tx.data.end(), rx.size = tx.data.size();
                ry.first = ty.data.begin(), ry.last = ty.data.end(), ry.size = ty.data.size();
                kar_rec_mul(rx, ry);
                sign = tx.sign != ty.sign ? +1 :  - 1;
            }
            z0.kar_rec_mul(x0, y0);
            add(z0);
            if(x.size >= n && y.size >= n){
                z2.kar_rec_mul(x1, y1);
                add(z2);
            }
            elemental_shift(n);
            add(z0);
            if(x.size >= n && y.size >= n){
                add(z2, n * 2);
            }
            normalize_data_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 unsigned_div(const integer &lhs, const integer &rhs){
            quo_rem qr;
            qr.quo.data.reserve(lhs.data.size());
            qr.rem.data.reserve(rhs.data.size());
            qr.rem.data.push_back(0);
            for(std::size_t i = lhs.data.size() * static_cast<std::size_t>(BitNum) - 1; i + 1 > 0; --i){
                if(!qr.rem.data.empty()){
                    qr.rem <<= 1;
                }else{
                    qr.rem.data.push_back(0);
                }
                qr.rem.data[0] |= (lhs.data[i / static_cast<std::size_t>(BitNum)] >> (i % static_cast<std::size_t>(BitNum))) & 1;
                if(unsigned_less_eq(rhs.data, qr.rem.data)){
                    qr.rem.unsigned_sub(rhs);
                    std::size_t t = i / static_cast<std::size_t>(BitNum);
                    if(qr.quo.data.size() < t + 1){ qr.quo.data.resize(t + 1); }
                    qr.quo.data[t] |= 1 << (i % BitNum);
                }
            }
            return qr;
        }

        static quo_rem div(const integer &lhs, const integer &rhs){
            quo_rem qr = unsigned_div(lhs, rhs);
            qr.quo.sign = qr.rem.sign = lhs.sign == rhs.sign ? +1 :  - 1;
            return qr;
        }

        static UInt ceil_pow2_single(UInt x){
            std::size_t n = 0, m = 0;
            while(x > 0){
                if((x & 1) > 0){
                    ++n;
                }
                ++m;
                x >>= 1;
            }
            if(n == 1){
                return static_cast<UInt>(1 << (m - 1));
            }
            return static_cast<UInt>(1 << m);
        }

        void elemental_shift(std::size_t n){
            data.insert(data.begin(), n, 0);
        }

        std::size_t bit_num() const{
            return bit_num(data.back());
        }

        static std::size_t bit_num(UInt x){
            std::size_t n = 0;
            while(x > 0){
                ++n;
                x >>= 1;
            }
            return n;
        }

        std::size_t finite_bit_shift_lsr(std::size_t n){
            std::size_t
                digit = n / static_cast<std::size_t>(BitNum),
                shift = n % static_cast<std::size_t>(BitNum),
                rev_shift = static_cast<std::size_t>(BitNum) - shift;
            UInt c = 0;
            for(std::size_t i = 0; i < data.size(); ++i){
                UInt x = data[i];
                data[i] = (x << shift) | c;
                if(rev_shift < BitNum){
                    c = x >> rev_shift;
                }else{
                    c = 0;
                }
            }
            if(c > 0){ data.push_back(c); }
            return digit;
        }

        void bit_shift_lsr(std::size_t n){
            elemental_shift(finite_bit_shift_lsr(n));
        }

        void bit_shift_rsl(std::size_t n){
            std::size_t
                digit = n / static_cast<std::size_t>(BitNum),
                shift = n % static_cast<std::size_t>(BitNum),
                rev_shift = BitNum - shift,
                size = data.size();

            if(digit > 0){
                for(std::size_t i = 0, length = size - digit; i < length; ++i){
                    data[i] = data[i + digit];
                }
                for(std::size_t i = 0; i < digit; ++i){
                    data[size - i - 1] = 0;
                }
            }

            UInt c = 0;
            for(std::size_t i = 0; i < data.size(); ++i){
                std::size_t j = size - i - 1;
                UInt x = data[j];
                data[j] = (x >> shift) | c;
                if(rev_shift < BitNum){
                    x = x << rev_shift;
                }else{
                    c = 0;
                }
            }
            normalize_data_size();
        }

        void normalize_data_size(){
            if(data.empty()){
                sign = +1;
                return;
            }
            std::size_t i = data.size() - 1;
            for(; ; ){
                if(i == 0 && data.front() == 0){
                    data.clear();
                    sign = +1;
                    return;
                }
                if(data[i] == 0){
                    --i;
                }else{
                    break;
                }
            }
            data.resize(i + 1);
        }

#define MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(assign_op, op, type) \
    integer &operator assign_op(const type &rhs){ *this = *this op integer(rhs); return *this; }

        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(+= , +, int);
        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(-= , -, int);
        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(*= , *, int);
        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(/= , /, int);
        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(%= , %, int);

        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(+= , +, unsigned int);
        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(-= , -, unsigned int);
        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(*= , *, unsigned int);
        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(/= , /, unsigned int);
        MULTI_PRECISION_ASSIGN_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(%= , %, unsigned int);
    };

#define MULTI_PRECISION_COMPARE_OPERATOR(op, type) \
    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
    bool operator op(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const type &rhs){ \
        return lhs op integer(rhs); \
    } \
    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
    bool operator op(const type &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){ \
        return integer(lhs) op rhs; \
    }

    MULTI_PRECISION_COMPARE_OPERATOR(==, int);
    //MULTI_PRECISION_COMPARE_OPERATOR(!=, int);
    MULTI_PRECISION_COMPARE_OPERATOR(<=, int);
    MULTI_PRECISION_COMPARE_OPERATOR(>=, int);
    
    MULTI_PRECISION_COMPARE_OPERATOR(==, unsigned int);
    //MULTI_PRECISION_COMPARE_OPERATOR(!=,  other.sign){ return true; }
            if(sign unsigned int);
    MULTI_PRECISION_COMPARE_OPERATOR(<=, unsigned int);
    MULTI_PRECISION_COMPARE_OPERATOR(>=, unsigned int);

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator +(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        integer<UInt, DoubleUInt, DoubleInt, BitNum> t(lhs);
        t.add(rhs);
        return t;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator -(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        integer<UInt, DoubleUInt, DoubleInt, BitNum> t(lhs);
        t.sub(rhs);
        return t;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator *(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        return integer<UInt, DoubleUInt, DoubleInt, BitNum>::mul(lhs, rhs);
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator /(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        return integer<UInt, DoubleUInt, DoubleInt, BitNum>::div(lhs, rhs).quo;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator %(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        return integer<UInt, DoubleUInt, DoubleInt, BitNum>::div(lhs, rhs).rem;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator <<(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, UInt n){
        integer<UInt, DoubleUInt, DoubleInt, BitNum> x(lhs);
        x.bit_shift_lsr(n);
        return x;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator >>(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, UInt n){
        integer<UInt, DoubleUInt, DoubleInt, BitNum> x(lhs);
        x.bit_shift_rsl(n);
        return x;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator +=(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        add(rhs);
        return *this;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator -=(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        sub(rhs);
        return *this;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator *=(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        *this = integer<UInt, DoubleUInt, DoubleInt, BitNum>::mul(lhs, rhs);
        return *this;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> &operator /=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        lhs = integer<UInt, DoubleUInt, DoubleInt, BitNum>::div(lhs, rhs).quo;
        return lhs;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator %=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        lhs = integer<UInt, DoubleUInt, DoubleInt, BitNum>::div(lhs, rhs).rem;
        return lhs;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator <<=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, std::size_t n){
        integer<UInt, DoubleUInt, DoubleInt, BitNum> x(lhs);
        x.bit_shift_lsr(n);
        lhs = std::move(x);
        return lhs;
    }

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator >>=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, std::size_t n){
        integer<UInt, DoubleUInt, DoubleInt, BitNum> x(lhs);
        x.bit_shift_rsl(n);
        lhs = std::move(x);
        return lhs;
    }

#define MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(op, type) \
    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator op(const integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const type &rhs){ return lhs op integer<UInt, DoubleUInt, DoubleInt, BitNum>(rhs); } \
    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum> \
    integer<UInt, DoubleUInt, DoubleInt, BitNum> operator op(const type &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){ return integer<UInt, DoubleUInt, DoubleInt, BitNum>(lhs) op rhs; }

    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(+, int);
    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(-, int);
    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(*, int);
    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(/ , int);
    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(%, int);

    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(+, unsigned int);
    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(-, unsigned int);
    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(*, unsigned int);
    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(/ , unsigned int);
    MULTI_PRECISION_OPERATOR_OVERLOAD_FOR_EXPLICIT_TYPE(%, unsigned int);

    template<class UInt, class DoubleUInt, class DoubleInt, DoubleUInt BitNum>
    std::ostream &operator <<(
        std::ostream &os,
        integer<UInt, DoubleUInt, DoubleInt, BitNum> rhs
    ){
        using integer = integer<UInt, DoubleUInt, DoubleInt, BitNum>;
        std::vector<std::string> rseq;
        integer lo(10);
        for(; rhs.data.size() > 0; rhs /= lo){
            integer temp = rhs % lo;
            if(temp != 0){
                rseq.push_back(std::to_string(temp.data[0]));
            }else{
                rseq.push_back(std::to_string(0));
            }
        }
        os << (rhs.sign > 0 ? "" : "-");
        if(!rseq.empty()){
            for(auto iter = rseq.rbegin(); iter != rseq.rend(); ++iter){
                os << *iter;
            }
        }else{
            os << "0";
        }
        return os;
    }
}

#endif // MULTI_PRECISION_INCLUDE_INTEGER_HPP