わりざんするアルゴリズム(その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表現でもそのようになります(二の補数表現でももちろん同様です)。補正前の商が正しい商より小さい奇数だったか大きい奇数だったかにより、その二つの表現のどちらになるかが決まります。

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

なんか最近割り算する回路を書いていろいろ知ったのでまとめていきます。

以下、割り算とは以下のものをさすことにします。

  1. 商と余りを両方求める
  2. 符号付きの場合、商はゼロへ丸め、余りは被除数と同じ符号とする

2は、数学的な定義とは異なります。C99で定義された除算の定義に一致します。

回路の説明というよりは、先週作った任意幅整数型lbl::uintN_t<N>lbl::intN_t<N>GitHub - lpha-z/lbl: header only fixed width integer library)を使ってC++アルゴリズムを書き起こしたものの紹介です。

回復法(restoring division)

小学校で習うような、筆算で求める方法です。

「回復」の名は、いったん除数を引いてみて負になってしまうようであれば除数を足すことで元の状態に"戻す"、というところからきているようです。 ソフトウェア的に実装するならともかく、ハードウェアで実現するのであれば、引く前の値を持っておけばよいだけなので、この名前はやや違和感があります。

#include "lbl/xintN_t.hpp"

template<std::size_t N>
constexpr lbl::uintN_t<N> safe_abs( lbl::intN_t<N> x ) {
    return x < 0 ? -x : x;
}

template<std::size_t N>
constexpr auto divider( const lbl::uintN_t<N> dividend, const lbl::uintN_t<N> divisor ) {
    lbl::uintN_t<N> remainder = 0;
    lbl::uintN_t<N> quotient = 0;
    for( int i = N-1; i >= 0; --i ) {
        quotient <<= 1;
        remainder <<= 1;
        remainder |= dividend & 1ull<<i ? 1 : 0;

        const lbl::uintN_t<N+1> test_sub = static_cast<lbl::uintN_t<N+1>>(remainder) - static_cast<lbl::uintN_t<N+1>>(divisor);

        const bool success = static_cast<lbl::intN_t<N+1>>(test_sub) >= 0;

        quotient |= success ? 1 : 0;
        remainder = success ? static_cast<lbl::uintN_t<N>>(test_sub) : remainder;
    }
    return std::pair {
        quotient,
        remainder
    };
}

template<std::size_t N>
constexpr auto divider( const lbl::intN_t<N> dividend, const lbl::intN_t<N> divisor ) {
    const bool quotient_sign = dividend < 0 ^ divisor < 0;
    const bool remainder_sign = dividend < 0;
    const auto [quotient, remainder] = divider_impl( safe_abs(dividend), safe_abs(divisor) );
    return std::pair {
        quotient_sign ? -static_cast<lbl::intN_t<N>>(quotient) : static_cast<lbl::intN_t<N>>(quotient),
        remainder_sign ? -static_cast<lbl::intN_t<N>>(remainder) : static_cast<lbl::intN_t<N>>(remainder)
    };
}


reminder |= dividend & 1ull<<i ? 1 : 0;の部分は、筆算で「次の桁をおろしてくる」部分に相当します。dividendは他の部分で使われないので破壊しても問題なく、dividend自体をシフトすることによって、ハードウェアは簡素化しますが、見やすさを優先しました。そのほか、最初の方はremainderの上位桁が使われていないことを利用すると、実はdividendと共存可能で、記憶素子をケチることができます。

要するに、引いてみても正だったら商として1を立て、そうでなければ0を立てる、というだけのアルゴリズムです。

引いてみて正かどうかを判定する部分は、オーバーフローしないように、結果を一桁多く求める必要があります(ボローが発生するかの確認が必要です)。reminderdivisor自体はN桁のまま計算してよいのですが、C++ではうまく書けなかったのでこのような形になっています。

for文で書いてある部分は、ハードウェアでは普通、複数サイクルかけて計算を行うことになるでしょう。 for文を回る回数はN回であり、各ループ内で一回加算器*1を使っています。つまり、1クロックに加算器が一回しか使えないような状況*2では、割り算を終わらせるのにNクロックかかるということになります。


符号付きの除算を行う場合、符号がついていると面倒なことになりますが、被除数・除数ともに絶対値を取って計算をした後、最後の段階で符号を正しく付け直すと簡単になります。

絶対値を取ったり負号をつけたりするのにも加算器が必要なことを考えると、素直に作るとN+2クロックかかるはずです。少し工夫すると、最終サイクルでの符号調整はループの中に組み込めるので、ループの最後の回では処理を変えることにより、余分な2クロックのうち1クロックは除去可能です。

また、ループの最初の回で立つ商が1であるパターンは限られているので、そのパターンを書き下せば、加算器を使わずに商を立てることができます。それと並行して絶対値を取る処理を行えば、余分なもう1クロックを削ることができます。そこまでする必要があるかは疑問ですが……。


回復法で除算を行う場合、加算器の結果がマルチプレクサの入力になりますが、典型的なFPGA (Field programmable gate array)に実装する場合、このことが大きな弱点になります。 典型的なFGPAの論理セルでは、組み合わせ回路を実現する LUT (Look-up table) の直後にフルアダーが配置してあります。マルチプレクサはLUTで実現、加算器はフルアダーで実現することになりますが、順序が逆なため二つの論理セルを使わないと実現できないことになります。

*1:加算器と言っていますが、専ら減算に使っていますね……。

*2:長大な桁数の加算器は(繰り上がりが連鎖するような、最悪な場合に)遅延が大きく、直列に複数回使うと1クロック内での最大遅延が大きくなりすぎ、クロック周波数が下がる原因となりえます。

半端なbit長の、固定幅の整数型

C99/C++11で追加された、固定幅の整数型(uint8_tint32_tのようなものたち)は特定のビット幅で計算したいときに大変便利です。何ビットの整数が要求されているかを明示することで移植性が高まります。

また、uint16_tは必ずオーバーフローせず、結果の上位ビットが無視されるだけの結果になるなど、ハードウェアに近いことをやりたいときにも便利です。しかし、この用途ではint13_tuint42_tが欲しくなることはないでしょうか。

大きめのbit幅の整数型を用い、適宜ビットマスクを掛けることで所望の動作をシミュレートすることはできますが、ユーザーコードで行うのは不便です。 当然そのようなライブラリが求められますが、少し探した感じではそのようなライブラリは見当たりませんでした。

そこで自作してみたというのが今日のお話です。

GitHubで公開してありますので、興味があれば使っていただけると嬉しいです。間違いなどの指摘も歓迎です。

github.com

以下、このライブラリを作るにあたって、移植性と効率性を両立するために工夫した点です。

下位bitを使ってシミュレートするのではなく、上位bitを使ってシミュレートする

単純な発想では、uint13_tをシミュレートするためには、uint16_tの下位13bitを用いる方法が考えられます。つまり、演算の後に常に0x1fffとのビットごと論理積を取りその結果を保持するという実装です。

この方法は、標準にある整数型に変換するときのコストが全くかからないというメリットがありますが、毎演算ごとにビットごと論理積が発生するというデメリットがあります(コンパイラがどのくらい見抜いてくれるかは確認していないですが……)。

そうではなく、あえてuint16_tの上位13bitを使うという方法を使っています。この方法を使うと、上位ビットが無視される状況が何の追加演算もなしに得られます。また、表現は単に定数倍(uint13_tuint16_tでシミュレートする場合は、3bit分なので8倍)されているだけであるので、大小比較は特に工夫せずにそのまま行って問題ありません。

代わりに、右シフトの場合は下位ビットが0になるように少しだけ補正が必要です。また、乗除算はビットマスクを掛ける処理の代わりにオフセットをずらす分の処理が少し増えます(ライブラリの実装では、左右のオペランドのどちらもオフセットを取り除き通常の表現にした後計算し、オフセットを戻すという処理を行っています。オフセットを取り除く処理とオフセットを戻す処理は打ち消す合うはずなので冗長ですが、コーナーケースが怖く、安全な実装に逃げました)。

ナイーブな実装と比べ、標準にある整数型との変換コストがかかる点だけは劣りますが、シフト演算一回分であり、頻繁な変換は普通は発生しないと考えられるため、問題ないものと考えました。

safe_abs

safe_absC++で安全に絶対値を求める - よーる で紹介した、未定義動作であるオーバーフローを起こさずに絶対値を取る関数です(詳細は以前の記事を参照ください)。

bit_cast

bit_castは移植性のある方法で、符号なし整数型を符号あり整数型(二の補数表現)に変換する(ビット表現を解釈しなおす)関数です。単にstatic_castを使う場合、変換後の値が負になる場合、処理系定義の動作となります。std::memcpyC++20で入る予定のstd::bit_castを用いても同じ動作が実現でき、そのような問題はなくなりますが、constexpr化できなくなります。

この関数では、uintN_tintN_tの変換において値が変わらない場合は移植性が保たれることを利用し、

  • 値が変わらない範囲(つまり、0x0000...から0x7fff...の範囲)にある場合はそのまま変換
  • 値が変わってしまう範囲(つまり、0x8000...から0xffff...の範囲)にある場合は、0x8000...を引いた値を求め、その値を変換(0x0000...から0x7fff...の範囲にあるので変換は移植性がある)したのち、0x8000...に相当する値(intN_tの最小値にあたる値)を加算する(この加算ではオーバーフローは発生しえない)

という方法で移植性を担保しています。最適化コンパイラを用いた場合、この処理は消失します。

算術シフト

符号あり整数型の右シフトの挙動は処理系定義となっています。x86にはsar (Shift arithmetic right) 命令が存在するため(他の多くのアーキテクチャにも算術シフト命令はあります)、多くの処理系では算術シフトとなるはずですが、やはり移植性のあるものではありません。 このライブラリを用いた場合、必ず算術シフトとなります。

なるべくオーバーヘッドを減らしたいため、少なくともx86のg++やclang++ではsar命令にコンパイルされるようなコードを書きたいところですが、うまくいきませんでした。

移植性よく、しかし既存のx86処理系にsar命令を推論させる方法は、試した限りでは存在しませんでした。コンパイラsar命令を出力するのは、

  1. 符号あり整数型を右シフトした場合
  2. 定数での除算の場合

に限られるようです。特に、以下のようなコードであっても最適化されてsar命令になることはありませんでした。

constexpr uint64_t sar1( uint64_t x ) noexcept {
    const uint64_t sign = x & 0x8000'0000'0000'0000;
    return sign | (x>>1);
}

1の場合を移植性よく使うことはできませんので、2の場合をうまく活用する方法を考えたいところですが、右算術シフトと除算では、被除数が負の時の丸め方向が異なります(右算術シフトは負の無限大への丸め、除算はゼロへの丸め)。そのため、正負で場合分けするコードが出現してしまいます。

この定数除算がsar命令を使ったコードに変換される動作は、コンパイラバックエンドで行われているようです。つまり、コンパイラフロントエンドでの情報(特に、未定義動作が起こる値が来ているかどうか)の情報が失われてしまっている点が厄介です。

定数除算がsar命令に変換されるのはほぼ定型文として変換されているようで、正の数が来ないはず、などの情報を活用できていません。コンパイラフロントエンドでの情報が残っていれば、定型文として変換した後にあり得ない分岐が取り除かれる最適化がかかるはずですが、失われているためにそれができていません。

おわりに

このライブラリはあまりビット効率のことを考えていないので、uintN_t<3>などとしても64bit使ってしまっています。8bitで十分なので、そういった点は改良の余地がありますが、あまり気にしないことにしました。

また、通常の整数と混ぜて使いたいためexplicitがほとんどつけていません。意図せぬ暗黙変換が発生してわけのわからないことになることが発生しそうですが、利便性を優先しました。この辺は好みが分かれるところですが、あまり新設ではないかもしれません。

再びになりますが、興味があればぜひ使ってみてください。間違いなどの指摘も歓迎です。

invSqrtをニュートン・ラフソン法で求める方法の性質

先週は平方根関数をライブラリとして実装する方法をひとつ紹介しました。その方法は、デフォルトの丸めモード(最近接丸め)では正しく動作しますが、方向丸めの場合にはうまくいかないことが判明しました。

その中で、逆数平方根ニュートン・ラフソン法で求める、invSqrt関数を実装しました。その実装は、丸め誤差に気を使っていない、ナイーブな実装になっています。それが原因で方向丸めの場合にうまくいかないのかを確かめるため、この関数の性質について調べます。

invSqrt関数の定義

以下で他の実装についても議論するので、invSqrt0とします。

float apploxInvSqrt( float x ) {
    static_assert( sizeof(float) == sizeof(uint32_t) );
    uint32_t i;
    std::memcpy( &i, &x, sizeof i );
    i = 0x5f3759df - ( i >> 1 );
    float f;
    std::memcpy( &f, &i, sizeof f );
    
    return f;
}

float invSqrt0( float x ) {
    // overflows if x <= 0x1.d1cc6p-129 
    assert( isnormal(x) && x > 0 );

    float f = apploxInvSqrt(x);
    f *= 1.5f - 0.5f * x * ( f * f );
    f *= 1.5f - 0.5f * x * ( f * f );
    f *= 1.5f - 0.5f * x * ( f * f );    
    return f;
}

逆数平方根ニュートン・ラフソン法で求める方法の性質

  • (無限精度で計算した場合、)以下の数式で求める場合、収束は単調になる。

{ \displaystyle
f_{i+1} = f_i \times (1.5 - 0.5xf_i^2)
}

invSqrt0実装の性質

  • invSqrt0の実装方法では、係数1.5f - 0.5f * x * (f*f)の計算を無限精度でやった後一度だけ丸めたとしても、返り値に0.5ulpを超える誤差が乗ることがある。

たとえば、x = 0x1.00127cp+0fの場合、初期予測は0x1.eebaa8p-1f、一回加速後は0x1.ff1084p-1f、二回加速後は0x1.ffecf4p-1fとなりますが、この時つぎの係数は正確に計算すると0x1.00004882ac8fcp+0程度になります。これを丸めると0x1.000048p+0fになるが、これを乗じると0x1.ffed84p-1fになってしまいます。 切り上げで丸めると0x1.00004ap+0fを得ることができますが、これを乗じた場合は0x1.ffed88p-1fになってしまいます。つまり、真の値0x1.ffed850031aa1dp-1を丸めた0x1.ffed86p-1になるような係数は、そもそも単精度浮動小数点数の中に存在しません。

これは、求めるべき値の仮数部が2に近く、精度が必要であるのに、係数の仮数部は1に近く、精度が一桁分足りなくなってしまっていることが原因です。

  • invSqrt0の実装方法では、複数回の丸めが原因で、返り値に1.9ulpを超えるような誤差が乗ることがある。

たとえば、x = 0x1.08fd12p+0fの場合、二回加速後は0x1.f73d9ap-1fになります。この時、真の係数は0x1.00001a813f05dp+0程度です。計算過程を追ってみると、f2*f2 == 0x1.eea19p-1fですが、これは0x1.eea190ff6052p-1を丸めたものであり、0.5ulp近く切り下げられてしまっています。さらに0.5f * x * (f2*f2) == 0x1.ffff94p-2fですが、これは0x1.ffff95fb03e8bp-2を丸めたものであり、再び0.5ulp近く切り下げられてしまいました。 次に、1.5f - 0x1.ffff94p-2fが計算され、結果は0x1.00001bp+0となるはずですが、偶数丸めが発動して0.5ulp切りあがります。複数回丸めてはいけないという典型例になっています。

全ての丸めが同じ方向に働いたため、誤差が大きくなってしまいました。 ちなみに、この例では、誤差は1.9485ulpに達しています。

  • 丸め誤差が原因で、何度加速しても収束しないことがある。
    f *= 1.5f - 0.5f * x * ( f * f );

の加速を複数回繰り返しても、誤差が原因で二つの値を振動することがあります。

丸め誤差を減らすために最後にstd::fmaを使うinvSqrt1実装の性質

float kasoku( float x, float y ) {
    float hh = 0.5f - (0.5f*x) * (y * y);
    float ret =  std::fma(y, hh, y);
    return ret;
}

float InvSqrt1( float x ) {
    // overflows if x <= 0x1.d1cc6p-129 
    assert( isnormal(x) && x > 0 );

    float f = apploxInvSqrt(x);
    f *= 1.5f - 0.5f * x * ( f * f );
    f *= 1.5f - 0.5f * x * ( f * f );
    f = kasoku( x, f );
    
    return f;
}

このコードは、逆数平方根の収束公式と丸め誤差conv_2nd_B関数を参考にしました。

  • 最終回の加速の計算でstd::fmaを一回使う方法を用いても、誤差は1ulpには収まらない。

  • 上記のコードを用いてもう一度加速すると却って誤差が増える。

  • 上記のコードを用いてさらにもう一度加速すると、誤差は1ulpに収まる*1

丸め誤差を取り除くために追加で二回も加速しないといけないようでは、もはや加速とは呼べないような気がします。

二つのfloatを用いて極力丸め誤差を減らすinvSqrt2実装の性質

float kasoku( float x, float y ) {
    float y2 = y*y;
    float yc = std::fma( y, y, -y2 );
    float  h = std::fma( -0.5f*x, y2, 0.5f );
    float hc = std::fma( -0.5f*x, yc, h );
    
    float ret = std::fma( y, hc, y );
    return ret;
}

float InvSqrt2( float x ) {
    // overflows if x <= 0x1.d1cc6p-129 
    assert( isnormal(x) && x > 0 );

    float f = apploxInvSqrt(x);
    f *= 1.5f - 0.5f * x * ( f * f );
    f *= 1.5f - 0.5f * x * ( f * f );
    f = kasoku( x, f );
    
    return f;
}

このコードは、逆数平方根の収束公式と丸め誤差conv_2nd_C関数を参考にした。二つのfloatを用いて48bit分の精度を確保する試みです。

  • この極力丸め誤差を減らした加速方法を三回目の加速に用いても、必ずしも真値を丸めた値になるとは限らない。ただし、誤差は0.5004ulpを超えない*2

  • それどころか、上記のコードを用いた場合でも、何度加速しても収束しない場合がある。

x = 0x1.13e070p+1fの場合、真の逆数平方根0x1.5cc0a9000000bp-1程度であり、48bit精度ではどちらに丸めてよいか判別がつきません。実際、0x1.5cc0a8p-1f0x1.5cc0aap-1の間で振動します。

 まとめ

調べてみましたが、あまり得るものがありませんでした……。最後に出てきた、0x1.13e070p+1fの真の逆数平方根0x1.5cc0a9000000bp-1であるという例はかなり凶悪ですね。 逆数平方根を正確に求めようとすると結構泥沼のようです。

mySqrt関数の改善

もっとも単純なinvSqrt0の方法でも誤差は2ulp以内であり、std::fesetround(FE_DOWNWARD);されている時に先週公開したmySqrt関数をこれを使っても、ごく一部の平方数に対する平方根をぴったり1ulp小さく算出してしまう以外の間違いは発生しないようです。

invSqrt1を用いるとそのような間違いの頻度はさらに下がります。invSqrt2を用いた場合、invSqrt0よりは減りますがinvSqrt1よりは増えるという謎現象が発生します。

*1:xが1付近の値であるときしか調べていません

*2:xが1付近の値であるときしか調べていません

IEEE浮動小数点数における平方根演算の精度に関する覚書

IEEE浮動小数点数における演算では、丸め誤差が不可避です。特に、複数回の演算を繰り返すと丸め誤差が積もっていき、正確な値と大きく離れた答えを得てしまうことがあります。しかし、次の演算については、(数学的に)正確な値を求めた後、一回だけの丸めが発生することが、IEEE標準で規定されています。

浮動小数点数演算のできるCPUであれば、普通、四則演算や積和演算を行う命令を持っています。 しかし、平方根を正確に計算する命令を持たない命令セットも存在します。

そのような場合、平方根関数はライブラリ実装となるわけですが、どのように実装すれば要求を満たせるのでしょうか?

C++std::sqrtは正確に計算しているのか?

結論

しています。

標準の丸めモード、つまり最近接丸め(ぴったり半分ならば偶数へ丸める)を用いた時、単精度浮動小数点数floatについては全ての値で正しく丸められていることを確認しました。

確認は、Wandboxのclang7.0.0で行いました。

その他の性質

std::fmaとの関係

また、任意の単精度浮動小数点数xに対して、次のことが成り立つようです。

std::sqrt(x)は、std::abs(std::fma(s,s,-x))を最も小さくする正のsである。

これは自明ではありません。丸め誤差が0.499999999ulpであるほうに丸めるより、-0.500000001ulpであるほうに丸める方がstd::abs(std::fma(s,s,-x))が小さくなる、という場合があり得るわけです。そのようなぎりぎりの位置に平方根が存在することがないことは証明できるのでしょうか。四則演算の場合はそんな感じの定理があったような……?

これは実験的に確かめたことであり、少し怪しいです。特に、xが非正規化数である場合はそのようなsが複数個あることがあります(s*s-xを無限精度で計算できる場合には、sが正かつfloatで表現可能、という条件のもとそのようなsは一意に定まるはずですが、std::fma(s,s,x)がほぼ0(絶対値が最小の正の非正規化数よりも小さく、丸めると0.0f-0.0fになる値)になり、精度が不足することが原因のようです)。

丸め誤差の最大値

丸め誤差が0.5ulp(丸める桁の0.5倍)にかなり近い状況は発生することがあります。 たとえばx = 0x1.000002p+0f1.0f + std::numeric_limits<float>::epsilon()とかstd::nextafter(1.0f,2.0f)などとも表せます。1+2^-23です)な場合にstd::sqrt(x) == 1.0fになりますが、その丸め誤差は0.499999985ulpになります。

平方根を取ってから二乗するとどうなるか

平方根を計算してから、その値を二乗すると元に戻らないことがあります。 たとえば、x = 0x1.ffd508p+0fに対して、std::sqrt(x) == 0x1.69fab4p+0f(真の値は大体0x1.69fab4fffe986p+0くらい)であり、std::sqrt(x)*std::sqrt(x) == 0x1.ffd506p+0fになりますが、左辺の結果を丸める前は0x1.ffd5052c0e90p+0であり、真の値から1.413951ulpもずれた値になっています。

しかし、このずれはほぼ√2 ulpに収まります。つまり、1.5ulpを超えることはないので、xstd::sqrt(x)*std::sqrt(x)の間の違いは最大でも1ulpに収まります。

二乗してから平方根を取るとどうなるか

なんと必ず元の値に戻ることが証明されているようです。

https://fixedpoint.jp/2016/03/18/sqrt-of-square-of-fp-number.html

もちろん、二乗した値が0になったり、INFINITYになった場合は元に戻ることはありません。

どういったアルゴリズムで計算できるのか

超越関数の場合、意地悪な引数を与えられると丸め誤差が0.5ulpに非常に近くなり、どちらに丸めてよいかを判断するために追加で計算する桁数が多くなる例があるようです。 平方根関数の場合、そのようなことはなく、仮数部の桁数から計算可能な定数桁を計算すれば丸める方向が決まるのでしょう(調べられませんでしたが、IEEE平方根関数の丸めは正確に行うと決められているのですから、何かしらの裏付けがあるのだと思います)。

単純な方法としては、整数上で十分な精度で計算してそれを浮動小数点数に変換するというものがあります。 しかし、この方法は明らかに遅く、実用に耐えるものとは言えません。 なんとか、浮動小数点演算をうまく活用して平方根関数を実現できないものでしょうか。

以下のような実装をしてみたところ、全ての浮動小数点数の値でstd::sqrtと値が一致しました。ただし、std::sqrtと違い、エラー報告をしない点は動作が違いますので注意してください。 負の数を入れると-nanが返ってくるんですね。

float apploxInvSqrt( float x ) {    
    static_assert( sizeof(float) == sizeof(uint32_t) );
    uint32_t i;
    std::memcpy( &i, &x, sizeof i );
    i = 0x5f3759df - ( i >> 1 );
    float f;
    std::memcpy( &f, &i, sizeof f );
    
    return f;
}

float invSqrt( float x ) {
    // overflows if x <= 0x1.d1cc6p-129 
    assert( isnormal(x) && x > 0 );

    float f = apploxInvSqrt(x);
    f *= 1.5f - 0.5f * x * ( f * f );
    f *= 1.5f - 0.5f * x * ( f * f );
    f *= 1.5f - 0.5f * x * ( f * f );    
    return f;
}

float mySqrt( float x ) {
    if( x == 0 ) { 
        // same sign
        return x;
    }
    if( x < 0 ) {
        return std::copysign( std::nanf(""), -1.0f );
    }
    if( std::isinf( x ) ) {
        return x;
    }
    if( x < std::numeric_limits<float>::min() ) {
        return std::ldexp( mySqrt( std::ldexp( x, 32 ) ), -16 );
    }
    const float q = invSqrt(x);
    const float t = std::fma( x*q, x*q, -x );
    const float X = x*q - std::ldexp( q*t, -1 );
    
    return X;
}

apploxInvSqrt関数は、Fast Inverse Square Rootと呼ばれるビット演算テクニックを使っています(C++ではビット列を読み替えるためにキャストやunionを使うのはstrict alias rulesに反しますので、std::memcpyを使います。C++20ではstd::bit_castというそのものずばりの関数が導入されるのでそのような間違いは減ってゆくでしょう)。

invSqrt関数では、apploxInvSqrtで得られる初期値をもとに、ニュートン・ラフソン法を用いて逆数平方根の精度を上げています。もともと4bitは精度があるため、3回加速すれば8bit, 16bit, 32bitと精度が上がり、24bit精度のfloatでは十分な精度になります。 逆数平方根の精度が重要な場合は、最後の加速は丁寧に行う必要がありますが(参考:逆数平方根の収束公式と丸め誤差)、今回は平方根を求めることができればいいため、細かいことは気にしない実装になっています。

mySqrt関数では、基本的にはmySqrt(x) = x * invSqrt(x)の式によって計算を行います(これは、低速な除算を避けるためです。平方根自体をニュートン・ラフソン法で求めるのではなく、逆数平方根ニュートン・ラフソン法で求めたのも除算を避けるためです。平方根ニュートン・ラフソン法で求めたい場合、除算が必要になってしまいます)。 これだけだと最後の一桁が正しくならないことがあるので、(一次の)ニュートン法で補正します。微分係数は1/2√xであり、二次収束が必要な場合は1/2(x*q)の計算が必要になりますが、一次収束で十分なのでq/2で代用できます。 なお、非正規化数の平方根を求める場合、逆数平方根を介すと大きくなりすぎて表現できなくなる危険性があります。そのため、非正規化数であった場合は232をかけた値の平方根を求め、それを216で割ることにより、この問題を回避しています。

標準でない丸めモードを使うとどうなるのか

最近接丸め(丸め誤差がぴったり0.5ulpになることは平方根演算では発生しないので、偶数丸めでも無限大への丸めでもどちらでもよい)の場合は上記の方法で求められることが分かりました。

方向丸め(正の無限大への丸め、負の無限大への丸め、零への丸め)の場合もうまくいくのでしょうか?

結論から言うと、上記の方法では正しく求められません。

平方根関数の丸めモードを指定する場合、以下のようなコードになるでしょう。

std::fesetround(FE_TOWARDZERO);
std::cout << std::sqrt(2.0f) << std::endl;

このように書いた場合、std::sqrt関数がライブラリ実装されており、その中で浮動小数点演算を行っている場合、その浮動小数点演算の丸めモードが意図せず変わってしまうことになります。

たとえば、零への丸めの場合、上記の方法ではx = 1.0fの時、q = invSqrt(x)についてq == 0x0.fffffep+0となってしまい、t = std::fma(x*q,x*q,-x)についてt == -1e-23になり、mySqrt(x) == 0x0.ffffffp+0となり、正しくありません。

tは正負両方の値を取りうるため、三つの丸めをすべて異なる方法で取り扱う必要が出てきます。いったいどうやったらよいのでしょう……?分岐してよいのなら、近似値を計算した後にstd::fmaで検算をして、ずれていれば1ulp増減させるという方法がありますが……。

やはりライブラリ実装をするのならば、整数上で計算するのが無難だということでしょうか。

しかしstd::sqrtは正しく計算するようです。しかも高速です。SSE命令のsqrtss命令を使っているからのようです。

そういえば先ほどの実装ではstd::fmaはちゃんと精度よく求まることを前提としていましたが、std::sqrtがライブラリ実装のアーキテクチャstd::fmaだけはハードウェア命令で実現されるなんてことはあるのかと思われるかもしれません。実際、Intel Itanium (IA-64)では、FMA命令はあるものの、FDIVやFSQRT命令は存在せず、そのようなアーキテクチャの一例になっています。

*1:ただし、-0に平方根演算を適用すると-0が返るらしい。

bit列のパターンマッチング分解を行う方法

先週の話題

bit列のマッチングを行う方法 - よーる

に引き続き、bit列操作の可読性を上げようという試みです。

パターンマッチを行い、特定のパターンに合致したと判明した時、残りの部分を(そのパターン特有の方法で)分解することは、特に関数型言語においては頻出の操作です。

bit列操作において頻出の操作であるかについてはやや疑問が残るところですが、使いたくなったので作ってみました。

まず、指定されたbitを抽出する、extract_bits関数を作ります。

#include <cstdint>
#include <utility>

constexpr uint64_t make_mask( const char* str, char C, uint64_t acc = 0 ) {
    return *str ? make_mask( str+1, C, acc<<1 | *str==C ) : acc;
}

template<uint64_t Mask, std::size_t... Indeces>
constexpr uint64_t extract_bits_impl( uint64_t bits, std::index_sequence<Indeces...> ) {
    uint64_t result = 0;
    int i = 0;

    ( (result |= Mask&1ull<<Indeces ? !!(bits&1ull<<Indeces)<<i++ : 0), ... );
    return result;
}

template<char C, class Str>
constexpr uint64_t extract_bits( Str pattern, uint64_t bits ) {
    return extract_bits_impl<make_mask( pattern(), C )>( bits, std::make_index_sequence<64>() );
}

#include <iostream>
int main() {
    uint64_t x;
    std::cin >> x;
    std::cout << extract_bits<'a'>( []{return"aaaabbbb";}, x ) << std::endl;
    std::cout << extract_bits<'b'>( []{return"aaaabbbb";}, x ) << std::endl;
    std::cout << extract_bits<'a'>( []{return"aaaaaaaabbbbbbbb";}, x ) << std::endl;
    std::cout << extract_bits<'b'>( []{return"aaaabbbbaaaabbbb";}, x ) << std::endl;
}

このコードは、実行時の効率にも注意を払っています。たとえば、for文で64回ループを作るのではなく、コンマ演算子の畳み込みを用いることで、ループアンローリングを強制します。 このようにすることで、コンパイラの定数伝搬最適化を助けることができます。

実際、このコードをclang++ -std=c++1z -O2 -S -masm=intel(バージョンは4.0.1、古い……)でコンパイルすると、以下のようなコードになります。

# extract_bits<'a'>( []{return"aaaabbbb";}, x )
mov   edx, dword ptr [rbp - 8]
shr   edx, 4
and   edx, 15

# extract_bits<'b'>( []{return"aaaabbbb";}, x )
mov   rdx, qword ptr [rbp - 8]
and   edx, 15

# extract_bits<'a'>( []{return"aaaaaaaabbbbbbbb";}, x )
mov   eax, dword ptr [rbp - 8]
movzx edx, ah  # NOREX

# extract_bits<'b'>( []{return"aaaabbbbaaaabbbb";}, x )
mov   rdx, qword ptr [rbp - 8]
mov   eax, edx
and   eax, 15
shr   rdx, 4
mov   ecx, edx
and   ecx, 16
or    rcx, rax
mov   eax, edx
and   eax, 32
or    rax, rcx
mov   ecx, edx
and   ecx, 64
or    rcx, rax
and   edx, 128
or    rdx, rcx

上三つのコードは最適コードになっています。最近のコンパイラはさすがの最適化性能です。最後のものは上位4bitの計算部分が最適とは言えなくなっていますが、抽出すべきbit部分が分離している時点で相当難しいことを要求しているため、これくらいは仕方がないでしょう(後半のor命令では不要に64bitレジスタを使っているところからもコンパイラが"見抜け"ていないことがうかがえます)。

パターンマッチを行うところまで実装すると、以下のようになります。

#include <cstdint>
#include <utility>
#include <array>

constexpr uint64_t make_mask( const char* str, char C, uint64_t acc = 0 ) {
    return *str ? make_mask( str+1, C, acc<<1 | *str==C ) : acc;
}

template<uint64_t Mask, std::size_t... Indeces>
constexpr uint64_t extract_bits_impl( uint64_t bits, std::index_sequence<Indeces...> ) {
    uint64_t result = 0;
    int i = 0;

    ( (result |= Mask&1ull<<Indeces ? !!(bits&1ull<<Indeces)<<i++ : 0), ... );
    return result;
}

template<char C, class Str>
constexpr uint64_t extract_bits( Str pattern, uint64_t bits ) {
    return extract_bits_impl<make_mask( pattern(), C )>( bits, std::make_index_sequence<64>() );
}



constexpr uint64_t make_mask( const char* str, uint64_t acc = 0 ) {
    return *str ? make_mask( str+1, acc<<1 | *str=='0' | *str=='1' ) : acc;
}

constexpr uint64_t make_bits( const char* str, uint64_t acc = 0 ) {
    return *str ? make_bits( str+1, acc<<1 | *str=='1' ) : acc;
}

template<class Str>
constexpr bool match( Str pattern, uint64_t test_bits ) {
    constexpr uint64_t mask = make_mask( pattern() );
    constexpr uint64_t bits = make_bits( pattern() );
    
    return (test_bits&mask) == bits;
}

template<char... Chars>
class Match {
    bool success;
    std::array<uint64_t, sizeof...(Chars)> decomp;
public:
    template<class Str>
    constexpr Match( Str pattern, uint64_t test_bits )
        : success( match( pattern, test_bits ) )
        , decomp { extract_bits<Chars>( pattern, test_bits )... }
        {}
        

    template<std::size_t N>
    constexpr auto get() const {
        if constexpr( N == 0 ) { return success; }
        else { return std::get<N-1>( decomp ); }
    }    
};

namespace std {
    template<char... Chars>
    class tuple_size<Match<Chars...>> : public std::integral_constant<std::size_t, sizeof...(Chars) + 1> {};
    
    template<std::size_t N, char... Chars>
    class tuple_element<N, Match<Chars...>> {
    public:
        using type = decltype( std::declval<Match<Chars...>>().template get<N>() );
    };
}

#include <iostream>
int main() {
    const uint64_t x = 0b00000000101000000000010110010011;
    auto [ADDI, rd, rs1, imm] = Match<'d','a','i'>( []{return"iiiiiiiiiiiiaaaaa000ddddd0010011";}, x );
    if( ADDI ) {
        std::cout << "ADDi x" << rd << ", x" << rs1 << ", " << imm << std::endl;
    }
}

[Wandbox]三へ( へ՞ਊ ՞)へ ハッハッ

ビット列のパターンマッチング分解なんて何に使うのかと思った方もいらっしゃったかもしれませんが、このように機械語のビット列をパースするのに使えます。 main関数では、RISC-VというISAのADDi x11, x0, 10という命令の機械語をパースしています(本来は符号拡張を実装しないといけないですが……)。

bit列のマッチングを行う方法

bit列操作のような低級なことをやっていると、bit列に対してマッチングを行いたいことがあります。 たとえば、以下のような例です。

t=0b11011010pattern=0b11xx1010にマッチするか?→マッチする。

このようなことを行うためには、pattern中に出現するx0に、それ以外を1に置き換えたmask=0b11001111pattern中に出現するx0に置き換えたbits=0b11001010を作成し、 (t&mask)==bitsであるかを確かめれば十分であり、一般的には最高速です(pattern=0b0000xxxx*1などの特殊な場合には、符号なし比較だけで判定できるため、常に最高速ではありません)。

さて、このようなmaskbitsを手作業で作成するのは発見しづらいミスのもとであり、ソースコードの意図も分かりづらくなります。 ソースコード中にpatternを文字列として記述できれば、そのような問題を解決できます。

しかし、そのような方法を用いた場合、実行時に処理されるため、低速になることが予想されます。 bit列操作のような低級なことをやっている場合、高速化のために行っていることが多く、いくらソースコードが見やすくなっても低速なのでは意味がないという状況も考えられます。

こういったときに便利なのがconstexprです。以下のように少しコードを書くだけで、上記の問題を解決することができます。以下のコードはC++17以降で動作します *2

#include <cstdint>
#include <cassert>

constexpr uint64_t make_mask( const char* str, uint64_t acc = 0 ) {
    assert( *str == '\0' || *str == '0' || *str == '1' || *str == 'x' );
    return *str ? make_mask( str+1, acc<<1 | *str!='x' ) : acc;
}

constexpr uint64_t make_bits( const char* str, uint64_t acc = 0 ) {
    assert( *str == '\0' || *str == '0' || *str == '1' || *str == 'x' );
    return *str ? make_bits( str+1, acc<<1 | *str=='1' ) : acc;
}


template<class Str>
constexpr bool match( Str pattern, uint64_t test_bits ) {
    constexpr uint64_t mask = make_mask( pattern() );
    constexpr uint64_t bits = make_bits( pattern() );
    
    return (test_bits&mask) == bits;
}

#include <iostream>
int main() {
    std::cout << match( []{return"xxxxxxxx1111xxxx";}, 0xacf9 ) << std::endl;
}

ところで、match関数の第一引数として、patternの文字列そのものではなく、ラムダ式を使っているのは、あまり普通の見た目ではないと言えるかもしれません。

マクロレス型安全printf - よーるでも使ったテクニックですが、文字列を返すラムダ式を使うことで、各々の文字列に異なる型を与えるといったことを行っています。

    constexpr const char* pattern = "xxxxxxxx1111xxxx";
    constexpr uint64_t mask = make_mask( pattern );
    constexpr uint64_t bits = make_bits( pattern );

のように自分でmaskbitsを計算してconstexpr変数に入れる場合は、ラムダ式を使った型付けテクニックを使わなくてもよいのですが、match関数のようにpatternも受け取る関数を作る場合はこのテクニックを使う必要があります。

パターンを文字列で書くのではなく、ユーザー定義リテラルで書く方法も考えられます。しかし、ユーザー定義リテラル(整数)は実際の整数リテラルとして有効な表現しか使えないという弱点があります。 今回の目的である可読性向上の目的には、今一歩使い勝手が悪いです。また、異なる型を与えないといけないという制約はユーザー定義リテラルを用いた場合でも満たさないといけないため、ユーザー定義リテラル(文字列)ではうまくいきません。ユーザー定義リテラル(文字列・テンプレート)があったらこんなラムダ式で囲むなんていう方法を取らなくてすむようになります。

*1:10/07訂正:もともとは0bxxxx1111と書いていましたが、誤りでした

*2:時代はC++17だというのに、いまだにループで回すみたいなコードではなく実質リターン文ひとつみたいなコードを書いてしまっています(別に悪いというわけではないはずですが)。あと-Wparenthesesに怒られますが、演算子の優先順位は把握している(し、スペーシングもそれを反映したものになっている)のになぁって思ってしまいます。まぁ優先順位を把握していないプログラマーに不親切だから、常識的にわかる優先順位以外はかっこでくくろうという意見も分からないではないですが。特にC++は暗黙変換で型検査を通ってしまう怖い言語なので、用心しようという姿勢は大切ですね(このコードでは暗黙変換を使っているわけですが)。