わりざんするアルゴリズム(その2.5)

三週間前の記事（わりざんするアルゴリズム(その2) - よーる）の割り算する回路・アルゴリズムを改良する話です。

わりざんするアルゴリズム(その1) - よーるとわりざんするアルゴリズム(その2) - よーるで紹介した方法は、商と剰余を求めるのに常に整数型のビット幅＋αサイクルかかっていました。つまり、9 / 7を計算するというように、実際に入っている値（の絶対値）が小さい場合でもそれだけのサイクル数がかかってしまうということです。

回復法の場合

回復法を使った除算法において、実際に立つ商が小さい場合、商の上位桁はずっと0が立ち続けることになります。筆算を行う場合、上位桁の0は普通書きません。つまり、上位桁に0を立て続ける部分は、実質的な計算が行われているわけではなく、無駄であるということになります。

実際にどの位に初めて商として1が立つのかは計算してみないとわからないですが、除数・被除数のleading zeros（最上位桁から連続する0の個数）を調べれば、初めて商として1が立つ桁か、その一つ上の桁を求めることはできます。これにより無駄な計算を極力省くことができます。

非回復法の場合

非回復法の場合、回復法のような単純な方法を使うことはできません。これは、商の表現が、各桁に $\overline{1}$ （上線付きの1であり、-1を表します。以下形の似ているTで代用します）と1しか使えない（つまり、0が使えない）ためです。

しかし、全く解決策がないわけではありません。商としてT111111111111TT1T == -27が立つ状況を考えてみます。立つ商が負でその絶対値が小さい場合、上位桁はT1111....という形になることが分かります（商が正でその絶対値が小さい時は、上位桁は1TTTT...という形になります）。この部分の計算はやや無駄であり、なんとか計算サイクル数を減らすことを考えます。

回復法での改良のように、途中から計算を始めることを考えます。もし上位桁に0が使えるのであれば、0000000000T1TT1T == -27のように表すことができます。したがって、三十二の位から商を立て始めればよいことになります。

その後、BSD表現から通常の二の補数表現に変換するところに二つの注意が必要です。

二の補数表現における符号拡張

回復法の場合は符号なしで計算していたため、計算していない上位桁は0で埋めればよかったのですが、非回復法の場合は符号付きで計算しているため、計算していない上位桁を何で埋めるのかが問題になります。

BSD表現では、最上位桁にTが立っていれば負、1が立っていれば正なので、それを利用して符号拡張用のマスクを作成します。

ビット演算テクニックによるBSD表現から二の補数表現への変換

もともとの非回復法では、quotient = (quotient << 1) | 1;というビット演算テクニックでBSD表現から二の補数表現へ変換していました。この時、最上位ビットが左シフトによって暗黙に追い出されること込みでこの変換は成り立っています。

途中から計算を始めた場合、商の長さが本来のビット幅より短くなっているため、未計算の上位桁に追い出されることになります。この追い出されたビットは不要であるため、適切にマスクを掛けて消すことが必要です。

実際、最上位桁がTの場合、それは商が負であることを示しているため、未計算部分は1で埋める必要があります。このとき、最上位桁のTは0で表されていますが、左シフトで追い出されてきたこの0も含めて1で埋める必要があります。また、最上位桁が1の場合、それは商が正であることを示しているため、未計算部分は0で埋める必要があります。このとき、左シフトで追い出されてきた最上位桁の1も含めて0で埋める必要があります。

この部分の手順を間違えると答えが完全におかしくなるので注意してください。

以上の二つを踏まえて改良したアルゴリズムをC++で記すと以下のようになります

改良されたアルゴリズムのC++コード

以下のコードは-std=c++17オプションをつけてコンパイルします。また、ライブラリとしてGitHub - lpha-z/lbl: header only fixed width integer libraryが必要です。

#include "lbl/xintN_t.hpp"

template<std::size_t N>
constexpr int countSignificantBits( lbl::intN_t<N> x ) {
    bool top_bit = x & 1ull<<(N-1);
    for( int i = N-2; i >= 0; --i ) {
        bool ith_bit = x & 1ull<<i;
        if( top_bit != ith_bit ) { return i+2; }
    }
    return 1;
}

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::intN_t<N> dividend, lbl::intN_t<N> divisor ) {
    lbl::intN_t<N> remainder = dividend < 0 ? -1 : 0;
    lbl:intN_t<N> quotient = 0; // any value is allowed

    int loop_count = countSignificantBits( dividend );
    for( int i = loop_count-1; i >= 0; --i ) {
        bool new_bit = dividend & 1ull<<i;
        if( (remainder < 0) != (divisor < 0) ) {
            quotient = (quotient<<1) | 0;
            remainder = ((remainder<<1) | new_bit) + divisor;
        } else {
            quotient = (quotient<<1) | 1;
            remainder = ((remainder<<1) | new_bit) - divisor;
        }
    }

    bool positive = quotient & 1ull<<(loop_count-1);
    quotient = (quotient << 1) | 1;
    if( positive ) {
       quotient &= (1ull << loop_count) - 1;
    } else {
        quotient |= ~((1ull << loop_count) - 1);
    }

    if( remainder == 0 ) {
        /* do nothing */
    } else if( remainder == divisor ) {
        assert( divisor < 0 );
        quotient += 1;
        remainder = 0;
    } else if( remainder == -divisor ) {
        assert( divisor > 0 );
        quotient -= 1;
        remainder = 0;
    } else if( (remainder < 0) != (dividend < 0) ) {
        if( (remainder < 0) != (divisor < 0) ) {
            quotient -= 1;
            remainder += divisor;
        } else {
            quotient += 1;
            remainder -= divisor;
        }
    }

    return std::pair {
        quotient,
        remainder
    };
} 

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::uintN_t<N> dividend, lbl::uintN_t<N> divisor ) {
    auto [quotient, remainder] = nonrestoring_divider<N+1>( static_cast<lbl::intN_t<N+1>>(dividend), static_cast<lbl::intN_t<N+1>>(divisor) );
    return std::pair {
        static_cast<lbl::uintN_t<N>>(quotient),
        static_cast<lbl::uintN_t<N>>(remainder)
    };
}

さらなる改良

先の方法では、商の最上位ビットは左シフトで追い出された後マスクを掛けられて消されます。しかし、符号拡張のための情報として利用されるため、商の最上位ビットが役に立っていないわけではありません。

そこで、商を立てるときにその処理をしてしまえば、後処理の部分で余計なことをする必要がなくなります。

具体的には、商の最上位ビットにTが立つ場合、それは商が負であることを意味しているため、符号拡張を見越してquotient = -1とします。商の最上位ビットに1が立つ場合、それは商が正であることを意味しているため、符号拡張を見越してquotient = 0とします。

それを実装したのが以下のC++コードになります。

#include "lbl/xintN_t.hpp"

template<std::size_t N>
constexpr int countSignificantBits( lbl::intN_t<N> x ) {
    bool top_bit = x & 1ull<<(N-1);
    for( int i = N-2; i >= 0; --i ) {
        bool ith_bit = x & 1ull<<i;
        if( top_bit != ith_bit ) { return i+2; }
    }
    return 1;
}

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::intN_t<N> dividend, lbl::intN_t<N> divisor ) {
    lbl::intN_t<N> remainder = dividend < 0 ? -1 : 0;
    lbl:intN_t<N> quotient = 0; // any value (or non-initiaziled) is allowed

    int loop_count = countSignificantBits( dividend );
    for( int i = loop_count-1; i >= 0; --i ) {
        bool new_bit = dividend & 1ull<<i;
        if( (remainder < 0) != (divisor < 0) ) {
            quotient = i == loop_count - 1 ? lbl::intN_t<N>(-1) : (quotient<<1) | 0;
            remainder = ((remainder<<1) | new_bit) + divisor;
        } else {
            quotient = i == loop_count - 1 ? lbl::intN_t<N>(0) : (quotient<<1) | 1;
            remainder = ((remainder<<1) | new_bit) - divisor;
        }
    }

    quotient = (quotient << 1) | 1;

    if( remainder == 0 ) {
        /* do nothing */
    } else if( remainder == divisor ) {
        assert( divisor < 0 );
        quotient += 1;
        remainder = 0;
    } else if( remainder == -divisor ) {
        assert( divisor > 0 );
        quotient -= 1;
        remainder = 0;
    } else if( (remainder < 0) != (dividend < 0) ) {
        if( (remainder < 0) != (divisor < 0) ) {
            quotient -= 1;
            remainder += divisor;
        } else {
            quotient += 1;
            remainder -= divisor;
        }
    }

    return std::pair {
        quotient,
        remainder
    };
} 

template<std::size_t N>
constexpr auto nonrestoring_divider( lbl::uintN_t<N> dividend, lbl::uintN_t<N> divisor ) {
    auto [quotient, remainder] = nonrestoring_divider<N+1>( static_cast<lbl::intN_t<N+1>>(dividend), static_cast<lbl::intN_t<N+1>>(divisor) );
    return std::pair {
        static_cast<lbl::uintN_t<N>>(quotient),
        static_cast<lbl::uintN_t<N>>(remainder)
    };
}