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

先週の記事(わりざんするアルゴリズム(その1) - よーる)に引き続き、割り算する回路・アルゴリズムの話です。

非回復法(non-restoring division)

回復法の場合、「計算してみて、その結果に応じて……」という、ある意味試行錯誤に近いアルゴリズムでした。 非回復法では、被除数(部分剰余)と除数の符号の関係のみによって商として何が立つかが決まり、そのような試行錯誤的な手順はありません。

さらに、以下の特徴があります。

  • 二の補数で表された符号付き整数を取り扱うことができる。
  • BSD表現という、特殊な表現で商を求める(最後に通常の二の補数表現に変換する)。
  • この手法だけでは商と余りを正しく求められないことがあるため、最後に補正をするステップが必要。

Binary signed digit (BSD)表現

コンピュータ上では、多くの場合二進法、つまり二を底とする位取り記数法を使って数値を表現しています。これは、各位が0か1で表され、その重みが2nとなっている記数法です。

BSD表現は、重みが2nである点は同じですが、各位の数が異なり、各位は-1か1で表されます。-1は \overline{1}のように、上線付きの1で表すこともあります。以降、見た目の似ている"T"で代用します。

まず、BSD表現は普通の意味での記数法ではない点に注意します。そもそも、奇数しか表せず、偶数が表せません。 しかし、表せるならその表現方法は一意です(冗長な記数法ではありません)。

一方、陽に符号のようなものを考えなくても負の数を表せる記数法である点は優れています(二の補数表現もそのような記数法の一つです)。 たとえばT111 (=-8+4+2+1=-1)など、最上位桁がTである数は負の数を表しています。

また、各位の表現として0が使えないため、「記されていない上位桁や小数点以下の部分は0」という暗黙の取り決めが適用できません。 二の補数表現の場合*1、「記されていない上位桁は、最上位桁と同じ数字が入る」という"符号拡張"を行うことで、任意の桁数に拡張することができます。 BSD表現の場合、以下のように最上位のTT11...1に、11TT...Tに、それぞれすることで任意の桁数に拡張することができます。

      T11T1 (=-16+8+4-2+1=-5)
T11111111T1 (=-1024+512+256+128+64+32+16+8+4-2+1=-5)

     1T1TTT (=32-16+8-4-2-1=17)
1TTTTTT1TTT (=1024-512-256-128-64-32-16+8-4-2-1=17)

C++コード

以下のコードは-std=c++17オプションをつけてコンパイルします。

#include "lbl/xintN_t.hpp"

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
    for( int i = N-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;
        }
    }

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

解説

非回復法の場合、回復法とは逆に符号付き整数を取り扱うのが基本であるため、符号なし整数の場合はN+1ビットの符号付き整数用割り算器を呼び出します。Nビットの符号付き整数は、N+1ビットの符号付き整数で表せるため、それを利用します。

非回復法の割り算器は大きく分けて二つのステージに分かれます。for文で書かれた一ビットずつ商を求める部分と、if-else文で書かれた商と余りを補正する部分です。

remainderの初期化は、以下でremainderの符号を確認する部分がありますので、それに備えて符号拡張された表現になるよう、0-1を入れておきます。 quotientの上位ビットは全く使われず、かつ以下の手順ですべてのビットが左シフトで追い出されるため何が入っていても構わないのですが、C++的には演算前に初期化しておかないと未定義動作ですので適当な値で初期化しました。ハードウェアに実装する場合は、変なゴミが入っていても問題なく動作するということです。

一ビットずつ商を求める部分では、余りと除数の符号を比較し、一致するか否かでそのサイクルでの動作が決まります。つまり、余りと除数の符号が同じであれば商として1を立てて余りから除数を引く、余りと除数の符号が異なれば商としてTを立てて余りに除数を足す、という動作です。Tはquotient変数の中では0として表現しておきます。

この動作は、その桁までで最近接丸めを行った商を求める動作になっています(回復法の場合、その桁までで切り捨てた商を求める動作になっているのと比較してください)。最近接丸めの場合、ぴったり真ん中だったときの動作が問題になります。それがどういう時に発生するかというと、その桁まででの余りが0になった場合に発生します。その場合でも何かの商を立てないといけないわけですが、0は二の補数表現における"符号"が"正"の側に属していますので、除数が正であれば1、除数が負であればTが立つことになります。

このように変な商を立てるのはまずいことになりそうですが、実際は安全です。ループ内では、 \rm -|divisor| \le remainder \lt |divisor|が不変条件として常に成立しています。

quotient = (quotient << 1) | 1;の部分は、BSD表現を通常の二の補数表現に変換する部分です。技巧的ですが、以下のように解釈することができます。BSD表現における、1になっている桁の寄与はquotient分です。また、Tになっている桁の寄与は、-(~quotient)分になります。つまり、BSD表現でのquotientは、二の補数表現ではquotient + (-(~quotient))に対応する値であるということが分かります。ここで、ビット演算テクニックを考えると、-(~q) == q+1であることが分かります。よって、quotient + (-(~quotient)) == quotient + quotient +1となります。これを単純化すると(quotient << 1) | 1になります。

最後に商と余りを補正する部分を見ていきます。この時点で、商は最近接奇数に丸めた値、余りはそれに対応する(商が奇数という制約の上での)絶対値最小剰余になっています。これをC99言語流の商と余りの制約に補正していきます。

  • 余りが0、つまり最終桁でぴったり割り切れた場合には補正は不要です。
  • 余りがぴったり除数と一致するというのは、本来商が偶数で割り切れるのが正しい時、かつ除数が負の場合に発生します。割り切れた次の位での動作では、商として0が立てられず、最近接丸めの動作から、除数が負の時にはTが立ちます。それ以下の桁では、常に余りと除数が一致し、商としてずっとTが立ちます(ループ内不変条件の等号が成立する状況です)。つまり、...T11111のような商が立っていることになります。この場合、1を足すことで正しい商と余りに戻します。
  • 余りがぴったり除数に負号をつけたものに一致するというのは、本来商が偶数で割り切れるのが正しい時、かつ除数が正の場合に発生します。その原理は上述のものと同一です。この場合1を引くことで正しい商と余りに戻します。
  • それ以外の場合、被除数と余りの符号が一致しなければ一致する側に戻します。あるいは、小数点以下の計算を行うことで、最近接偶数を求めているという解釈もできます(割り切れず、かつ商として最近接奇数が不適切ならば、最近接偶数は一意でそれが適切な商となります)。BSD表現では偶数は表せませんが、小数点以下の計算を無限に行えば、循環小数として表現可能です。偶数であれば小数点以下はすべて同じ数字が表れます*2。そのため、小数点以下第一位の桁に立つ商を求めればよいことになります。そのため、(remainder < 0) != (divisor < 0)という、商を求めるforループで使った条件式と同じものが再出現します。

特性

回復法と同様、一ビットずつ商を求めている特性から、ビット長Nに対して少なくともNクロックサイクルはかかる除算器になります。ただし、補正の部分を商を求めているサイクルと同一サイクルで行うのは無理があるので、普通はN+1クロックサイクルかかるはずです。また、符号なし整数の場合、N+1ビットの符号付整数としてこの除算器に入力するため、N+2クロックサイクルかかることになります。

FPGAに実装する場合、ループ部分の回路は典型的なFPGAの論理セルと非常に相性が良いです。典型的なFPGAの論理セルは、LUT→FA→FFという構成になっています。ここで、LUTはLook-up-tableの略で、組み合わせ回路を実現する素子、FAはFull adderの略で加算器、FFはFlip-flopの略で記憶素子です。if( (remainder) < 0) != (divisor < 0) )の部分がLUTに対応、その後の加算がFFに対応、その後に実行される結果を代入する部分がFFに対応、となっているためです。

*1:二の補数表現では各位の表現として0が使えますが、最上位桁の重みだけ -2^{N-1}と一様でないため、任意の桁数に拡張する際にはやはり暗黙のルールが適用できず、特別な操作が必要です。

*2:実数を小数点以下が無限に続く位取り記数法で表現すると1.0000...=0.9999...のように二通りの表現がありえますが、BSD表現でもそのようになります(二の補数表現でももちろん同様です)。補正前の商が正しい商より小さい奇数だったか大きい奇数だったかにより、その二つの表現のどちらになるかが決まります。