integer (multi-precision).

integer.hpp

#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 <random>
#include <type_traits>
#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{
        static_assert(sizeof(UInt) * 2 <= sizeof(DoubleUInt), "sizeof(UInt) * 2 <= sizeof(DoubleUInt)");
        static_assert(std::is_signed<DoubleInt>::value, "std::is_signed<DoubleInt>::value");

    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(data.empty() && other.data.empty()){ return false; }
            if(data.size() > 0 && other.data.empty()){ return false; }
            if(data.empty() && other.data.size() > 0){ return true; }
            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 *= -1;
                }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();
        }

        void unsigned_single_sub(UInt c, std::size_t i = 0){
            for(; c > 0; ++i){
                UInt t = data[i] - c;
                if(data[i] < c){
                    c = 1;
                }else{
                    c = 0;
                }
                data[i] = t;
            }
        }

        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{
            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;
        };

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

        void kar_mul_range(const kar_range &lhs, const kar_range &rhs){
            if(lhs.size == 0 || rhs.size == 0){
                data.clear();
                sign = +1;
                return;
            }
            data.resize(lhs.size + rhs.size);
            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(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last) const{
            return kar_less_impl<false>(rhs_first, rhs_last);
        }

        bool kar_less_eq(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last) const{
            return kar_less_impl<true>(rhs_first, rhs_last);
        }

        template<bool Final>
        bool kar_less_impl(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last) const{
            return kar_less_impl<Final>(data.begin(), data.end(), rhs_first, rhs_last);
        }

        template<bool Final>
        static bool kar_less_impl(
            typename data_type::const_iterator lhs_first, typename data_type::const_iterator lhs_last,
            typename data_type::const_iterator rhs_first, typename data_type::const_iterator rhs_last
        ){
            std::size_t lhs_n = std::distance(lhs_first, lhs_last), rhs_n = std::distance(rhs_first, rhs_last);
            if(lhs_n == 0 && rhs_n == 0){
                return false;
            }else if(lhs_n == 0 && rhs_n > 0){
                return true;
            }else if(rhs_n == 0 && lhs_n > 0){
                return false;
            }
            typename data_type::const_iterator 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 Final;
        }

        void kar_add(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last, std::size_t n = 0){
            std::size_t rhs_size = std::distance(rhs_first, rhs_last);
            if(rhs_size == 0){
                return;
            }
            rhs_size += n;
            if(data.size() < rhs_size){
                data.resize(rhs_size);
            }
            typename data_type::iterator operand_iter = data.begin() + n;
            DoubleUInt c = 0;
            typename data_type::const_iterator iter = rhs_first, end = rhs_last;
            if(iter != end){
                --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;
            }
            if(operand_iter != data.end()){
                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_signed_add(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last){
            if(sign < 0){
                if(kar_less_eq(rhs_first, rhs_last)){
                    integer s(std::move(*this));
                    data.assign(rhs_first, rhs_last);
                    unsigned_sub(s);
                    sign = +1;
                }else{
                    unsigned_sub(rhs_first, rhs_last);
                }
            }else{
                kar_add(rhs_first, rhs_last);
            }
        }

        void kar_sub(const typename data_type::const_iterator &rhs_first, const typename data_type::const_iterator &rhs_last){
            if(kar_less_eq(rhs_first, rhs_last)){
                integer s(std::move(*this));
                data.assign(rhs_first, rhs_last);
                unsigned_sub(s);
                sign *= -1;
            }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_rec_mul(kar_range(tx.data.begin(), tx.data.end()), kar_range(ty.data.begin(), ty.data.end()));
                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){
            return unsigned_div(lhs.data.begin(), lhs.data.end(), rhs.data.begin(), rhs.data.end());
        }

        static quo_rem unsigned_div(
            typename data_type::const_iterator l_first,
            typename data_type::const_iterator l_last,
            typename data_type::const_iterator r_first,
            typename data_type::const_iterator r_last
        ){
            quo_rem qr;
            std::size_t ln = std::distance(l_first, l_last), rn = std::distance(r_first, r_last);
            qr.quo.data.reserve(ln);
            qr.rem.data.reserve(rn);
            for(std::size_t i = ln * 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] |= (l_first[i / static_cast<std::size_t>(BitNum)] >> (i % static_cast<std::size_t>(BitNum))) & 1;
                if(kar_less_impl<true>(r_first, r_last, qr.rem.data.begin(), qr.rem.data.end())){
                    qr.rem.unsigned_sub(r_first, r_last);
                    std::size_t t = i / static_cast<std::size_t>(BitNum) + 1;
                    if(qr.quo.data.size() < t){
                        qr.quo.data.resize(t);
                    }
                    qr.quo.data[t - 1] |= 1 << (i % static_cast<std::size_t>(BitNum));
                }
            }
            return qr;
        }

        static quo_rem kar_rec_div(kar_range lhs, kar_range rhs){
            quo_rem qr;
            std::size_t n = (rhs.size + (rhs.size % 2)) / 2;

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

            {
                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;
                    }
                }

                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);
                    }
                }

                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);
                    }
                }

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

#if 0
        static quo_rem kar_large_div(kar_range lhs, kar_range rhs){
            quo_rem qr;
            integer &quo(qr.quo), &rem(qr.rem);
            quo.data.reserve(lhs.size);
            rem.data.assign(lhs.first, lhs.last);

            using iter_type = typename data_type::const_iterator;
            auto compare = [](iter_type rem_first, iter_type rem_last, iter_type rhs_iter, iter_type rhs_last){
                for(auto iter = rem_last - 1, jter = rhs_last - 1; ; --iter, --jter){
                    if(*iter > *jter){
                        return true;
                    }else if(*iter == *jter){
                        if(iter == rem_last){
                            break;
                        }
                    }else{
                        return false;
                    }
                }
                return true;
            };

            std::size_t i = rem.data.size() - 1;
            for(; rem.data.size() > rhs.size || (rem.data.size() == rhs.size && compare(rem.data.begin(), rem.data.end(), rhs.first, rhs.last)); --i){
                kar_range rem_range(rem.data.begin() + i, rem.data.end());
                quo_rem st = kar_rec_div(rem_range, rhs);
                quo.kar_add(st.quo.data.begin(), st.quo.data.end(), i);
                rem.data.resize(i + st.rem.data.size());
                std::copy(st.rem.data.begin(), st.rem.data.end(), rem.data.begin() + i);
            }

            return qr;
        }
#endif

        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 quo_rem modern_div(const integer &lhs, const integer &rhs){
            quo_rem qr = kar_rec_div(kar_range(lhs.data.begin(), lhs.data.end()), kar_range(rhs.data.begin(), rhs.data.end()));
            qr.quo.sign = qr.rem.sign = lhs.sign * rhs.sign;
            return qr;
        }

        // Miller Rabin's Primality Test.
        static bool miller_rabins_test(const integer &n, std::size_t s){
            static std::mt19937 g(0xABABABABu);
            integer a, m = n - 2;
            for(std::size_t i = 0; i < s; ++i){
                if(witness(m.random(g) + 1, n)){
                    return false;
                }
            }
            return true;
        }

        static bool witness(const integer &a, const integer &n){
            integer b = n - 1, d = 1, x, t;
            for(int i = static_cast<int>(b.bit_num()) - 1; i >= 0; --i){
                x = d;
                d = (d * d) % n;
                if(d == 1 && x != 1 && x != b){
                    return true;
                }
                t = b >> static_cast<UInt>(i);
                if(t > 0 && t.data[0] % 2 == 1){
                    d = (d * a) % n;
                }
            }
            return d != 1;
        }

        // Random Number
        template<class Generator>
        integer random(Generator &g) const{
            std::size_t s = g() % data.size() + 1;
            integer x;
            x.data.resize(s);
            for(std::size_t i = 0; i < s - 1; ++i){
                x.data[i] = g() & base_mask;
            }
            if(s == data.size()){
                x.data.back() = g() % data.back();
            }
            x.normalize_data_size();
            return x;
        }

        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()) + BitNum * (data.size() - 1);
        }

        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);
        }

        static typename data_type::const_iterator normalize_range(typename data_type::const_iterator first, typename data_type::const_iterator last){
            do{
                --last;
                if(last == first){
                    return last;
                }
            }while(*last == 0);
            return last + 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<UInt, DoubleUInt, DoubleInt, BitNum>(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<UInt, DoubleUInt, DoubleInt, BitNum>(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(<, int);
    MULTI_PRECISION_COMPARE_OPERATOR(>, int);

    MULTI_PRECISION_COMPARE_OPERATOR(==, unsigned int);
    //MULTI_PRECISION_COMPARE_OPERATOR(!=, unsigned int);
    MULTI_PRECISION_COMPARE_OPERATOR(<=, unsigned int);
    MULTI_PRECISION_COMPARE_OPERATOR(>=, 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 +=(integer<UInt, DoubleUInt, DoubleInt, BitNum> &lhs, const integer<UInt, DoubleUInt, DoubleInt, BitNum> &rhs){
        lhs.add(rhs);
        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.sub(rhs);
        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>::mul(lhs, rhs);
        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).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){
        lhs.bit_shift_lsr(n);
        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){
        lhs.bit_shift_rsl(n);
        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