分岐を増やすと速くなる不思議な現象

マイクロベンチマークを書いていたのですが、表題のようなことが起こったのでメモです。

環境

  • g++ (GCC) 12.1.0

問題のソースコード

a.cpp↓

#include <cmath>
#include <bit>
#include <tuple>

std::tuple<unsigned long, double, double> f(double w) {
    unsigned long n = std::bit_cast<unsigned long>(w) & 3;
    double x = std::fma(1.0, 1.0, w);
    double y = std::fma(1.0, 1.0, x);
    double z = std::fma(1.0, 1.0, y);
    return {n, y, z};
}

double g(double y, unsigned long n) {
    double ret = 1 + y * (2 + y * (3 + y * (4 + y * (5 + y))));
    return (n & 2) ? -ret : ret;
}

double bench(double w) {
    auto [n, y, z] = f(w);
    return g(y, n);
}

int dummy;

int main(int argc, char* argv[]) {
    for(size_t i=0; i < 100000000; ++i) {
        double w = std::bit_cast<double>((unsigned long)rand() << 24) * 1e300 * 1e10;
        double v = bench(w);
        if (dummy & 1) { dummy += v; }
    }
}

b.cpp↓

#include <cmath>
#include <bit>
#include <tuple>

std::tuple<unsigned long, double, double> f(double w) {
    unsigned long n = std::bit_cast<unsigned long>(w) & 3;
    double x = std::fma(1.0, 1.0, w);
    double y = std::fma(1.0, 1.0, x);
    double z = std::fma(1.0, 1.0, y);
    return {n, y, z};
}

double g(double y, unsigned long n) {
    double ret;
    if (n == 0)      ret = 1 + y * (2 + y * (3 + y * (4 + y * (5 + y))));
    else if (n == 1) ret = 1 + y * (2 + y * (3 + y * (4 + y * (5 + y))));
    else             ret = 1 + y * (2 + y * (3 + y * (4 + y * (5 + y))));
    return (n & 2) ? -ret : ret;
}

double bench(double w) {
    auto [n, y, z] = f(w);
    return g(y, n);
}

int dummy;

int main(int argc, char* argv[]) {
    for(size_t i=0; i < 100000000; ++i) {
        double w = std::bit_cast<double>((unsigned long)rand() << 24) * 1e300 * 1e10;
        double v = bench(w);
        if (dummy & 1) { dummy += v; }
    }
}

違いは、g関数の中身だけです。

$ diff a.cpp b.cpp
15c15,18
<     double ret = 1 + y * (2 + y * (3 + y * (4 + y * (5 + y))));
---
>     double ret;
>     if (n == 0)      ret = 1 + y * (2 + y * (3 + y * (4 + y * (5 + y))));
>     else if (n == 1) ret = 1 + y * (2 + y * (3 + y * (4 + y * (5 + y))));
>     else             ret = 1 + y * (2 + y * (3 + y * (4 + y * (5 + y))));

さて、a.cppとb.cppはどちらが速いでしょうか。

普通に考えると、分岐の少ないa.cppの方が速いように思われます*1。 少なくとも、a.cppとb.cppの速度が同じことはあっても、b.cppの方が速いということはありえなさそうに思われます。 しかし、実際に測ってみると、b.cppの方が圧倒的に速いです(rand()を呼ぶだけのコードとほとんど同じ速さです)。 これはなぜでしょうか。

両者の違い

両者の違いは、インライン展開されるかどうかにあります。 a.cppはbench関数が小さいため、main関数の中にインライン展開されます。 一方、b.cppはbench関数がやや大きいため、main関数の中にはインライン展開されません。

速くなる原因

速くなる原因は、bench関数の結果であるvを使っていないことにあります。

a.cppでは、main関数にインライン展開されたbench関数の中身が、dummy & 1の検査の前に実行されるようなコードにコンパイルされるようです。 一方、b.cppでは、dummy & 1の検査がbench関数の呼び出しより前に行われるコード、つまりbench関数の中身が実行されないようなコードにコンパイルされるようです。

つまり、a.cppはvを使っていないにもかかわらず、運よくbench関数の中身が実行されているだけです。 b.cppはvを使っていないので順当にbench関数を実行しない最適化が行えています。

まとめ

ベンチマーク関数の結果を真には使わないようなマイクロベンチマークを書いた場合、最適化でベンチマークコードが消されないかはそのベンチマークコードの大きさによることがあります。 今回の例では、インライン展開されないような大きなコードの場合のみ働く最適化がコードを消したことで、見かけ上速度が向上していました。

*1:b.cppはthen部とelse部で仕事が同じなので分岐がなくなりそうにも見えますが、g++の場合、n==0やn==1が成立する場合はn&2が成立しないことを利用した最適化が優先され、then部とelse部のマージは行われません。

間違ったポラードのロー法の使い方

概要

アルゴ式が公開している素因数分解するアルゴリズム(と称しているもの)は、特定の入力に対して正しく動作しません(無限ループします)。 ポラードのロー法が失敗した場合、疑似乱数生成器のシードのみを変更するのではなく、疑似乱数生成器自体を変更してリトライする必要があります。

アルゴ式が公開しているコードとうまくいかない入力

アルゴ式が公開しているコード

アルゴ式が公開している素因数分解する(と称している)コードでは、以下のようにポラードのロー法を使っています。

long long pollard(long long N) {
    if (N % 2 == 0) return 2;
    if (is_prime(N)) return N;

    auto f = [&](long long x) -> long long {
        return (__int128_t(x) * x + 1) % N;
    };
    long long step = 0;
    while (true) {
        ++step;
        long long x = step, y = f(x);
        while (true) {
            long long p = gcd(y - x + N, N);
            if (p == 0 || p == N) break;
            if (p != 1) return p;
            x = f(x);
            y = f(f(y));
        }
    }
}

この方法では、ポラードのロー法により非自明な因数を見つけられなかった場合(疑似乱数が真に循環してしまった場合)に、++stepでリトライしています。 ここで、stepの変更によりx、つまり疑似乱数生成器のシードを変更しています。

しかし、これは問題があります。 なぜなら、疑似乱数生成器が弱い(循環周期が短い)場合、シード値を変更してもその状況が変わらないからです。

うまくいかない入力

上記pollard関数に124376107291を入力すると、無限ループします。 この数は、352523×352817で表される半素数素数二つの積)です。

この数の素因数分解に失敗するのは、疑似乱数生成器  f : x \mapsto x^2+1\!\!\!\mod\!124376107291 が821という非常に短い周期を持つことが関係しています。 具体的には、 \forall x\in \mathbb{Z}, f^{653}(x) = f^{1474}(x) が成り立ちます。 この短い周期により、mod 124376107291で循環する前にmod 352523での循環やmod 352817での循環を発見することができなくなっています。

どうしたらよいのか

ポラードのロー法が失敗した場合、疑似乱数生成器のシードを変えるのではなく、疑似乱数生成器自体を変えるほうが良いです。 具体的には、以下のように変更します。

      if (N % 2 == 0) return 2;
      if (is_prime(N)) return N;
  
+     long long step = 0;
      auto f = [&](long long x) -> long long {
-         return (__int128_t(x) * x + 1) % N;
+         return (__int128_t(x) * x + step) % N;
      };
-     long long step = 0;
      while (true) {
          ++step;

この変更を加えた場合、1013以下の数の非自明な因数を一つ見つけるのに必要なgcdの回数は最悪でも9000回になります。 1018以下の場合、最悪でもおおよそ12万回程度のgcd呼び出しで非自明な因数を一つ見つけることができるようです(私が見つけられた数の中で最も多くgcdを呼び出したのは、814483663644399613を入力したときで、117480回でした。2.1GHzで動く計算機で、25ms程度かかります)。

gcdの中で再帰的に呼び出されるgcdの呼び出し回数は含んでいません。

その他

乱数生成器のパラメータを変えるだけでは、問題を解決することはできません。

乱数生成器として  f : x \mapsto x^2 +2\!\!\!\mod\!N を使う場合、273772559の素因数分解に失敗します(無限ループします)。 乱数生成器として  f : x \mapsto x^2 +4\!\!\!\mod\!N を使う場合、2059の素因数分解に失敗します(無限ループします)。 乱数生成器として  f : x \mapsto x^2 +5\!\!\!\mod\!N を使う場合、385515865499の素因数分解に失敗します(無限ループします)。

乱数生成器として  f : x \mapsto x^2 +3\!\!\!\mod\!N を使う場合、無限ループする入力を見つけることはできませんでしたが、7816550168663の素因数分解には一分以上時間がかかります。

乱数生成器として  f : x \mapsto x^2\!\!\mod\!N を使う場合、無限ループすることはありません。 ただし、実質的に試し割り法に退化しているので、まったく実用的ではありません。 例えば7482809861の素因数分解に一分以上時間がかかります。

SRT法で平方根を計算する

先々週の記事(SRT除算について - よーる)で取り扱ったように、SRT法は除算をするアルゴリズムですが、少しだけ変形すると平方根を求めることにも使えます。

まずは、長除法(割り算を筆算で計算する手順)と開平法(平方根を筆算で計算する手順)の類似性を見ていきます。 以下では、基数を rとします。手で行う開平法の基数は r = 10です*1

長除法では、以下の手順を繰り返します。

  • 部分剰余 R_iと除数 Dから、部分商 Q_iを何らかの方法で求める。
  • 部分剰余を R_{i+1} = (R_i - Q_iD)rで更新する
    • 被除数 Nと暫定商 S_i = \sum_{k=0}^i Q_kr^{-k}の間に S_iD + R_{i+1}r^{-i-1} = Nが成り立つことを維持しながらアルゴリズムが進行する

開平法は、以下の手順を繰り返します。

  • 部分剰余 R_iと暫定商 S_{i-1} = \sum_{k=0}^{i-1} Q_kr^{-k}から、部分商 Q_iを何らかの方法で求める。
  • 部分剰余を R_{i+1} = (R_i - 2Q_iS_{i-1}-(Q_i)^2r^{-i})rで更新する。
    • 与数 Nと暫定商 S_i = \sum_{k=0}^i Q_kr^{-k}の間に (S_i)^2 + R_{i+1}r^{-i-1} = Nが成り立つことを維持しながらアルゴリズムが進行する

部分剰余の更新部分を見てみると、開平法では (Q_i)^2r^{-i}という余計な項がついているものの、この項を無視すれば長除法における Dと開平法における 2S_{i-1}が対応していることが分かります。 この対応は部分商 Q_iを求めるときの手順にも表れています。 つまり開平法は、その手順中で除数 Dを使う代わりに暫定商 S_{i-1}を使う(毎回“除数”が変化する)という点を除いて、ほとんど長除法と同じです。 よってSRTアルゴリズムを用いることで、平方根を求めることができます。

さて、 (Q_i)^2r^{-i}という余計な項を無視したのが気になります。 実はこの余計な項は、PDプロット上の境界線を少し上にずらす効果があります。 以下では、PDプロット上の境界線の表式を具体的に計算していきます。

PDプロット上の境界線の表式

SRTアルゴリズムを無限に続けていくことを考えます。 この時、暫定商が与数の平方根に限りなく近づかないといけません。 これは、残差 w_j = N - (S_j)^2が限りなく0に近づく、つまり \lim\limits_{j\rightarrow\infty}w_j = 0が成り立つ必要がある、と言い換えられます。

これが成り立つために、途中での残差 w_jはどのような条件を満たしている必要があるでしょうか。  w_\infty w_jとそれ以外に分解してみましょう。  S_\infty = S_j + \sum\limits_{k=j+1}^{\infty}Q_kr^{-k}なので  w_\infty = N - (S_\infty)^2 = N - \left({}(S_j)^2 + 2S_j\sum\limits_{k=j+1}^{\infty}Q_kr^{-k} + \left(\sum\limits_{k=j+1}^{\infty}Q_kr^{-k}\right)^2\right)となり、  w_j = N - (S_j)^2を使ってまとめると w_\infty = w_j - 2S_j\sum\limits_{k=j+1}^{\infty}Q_kr^{-k} - \left(\sum\limits_{k=j+1}^{\infty}Q_kr^{-k}\right)^2のようになります。  w_\inftyは0になるべきですから、 w_j = 2S_j\sum\limits_{k=j+1}^{\infty}Q_kr^{-k} + \left(\sum\limits_{k=j+1}^{\infty}Q_kr^{-k}\right)^2が成り立っている必要があります。 つまり考えるべきなのは、 Q_{j+1}, Q_{j+2}, \dotsをうまく選べば左辺と一致させられるか否かです。 右辺を最も大きくできるのは、各桁として許される数字の中で最も大きい数字を選び続けた時です。 右辺を最も小さくできるのは、各桁として許される数字の中で最も小さい数字を選び続けた時です。 その間の数であればどんな数でも表現できる(=記数法である)ことを確認する必要はありますが、常識的にそれは成り立っているとして、以下では上限と下限を計算します。

基数 r = 4で、各桁として許される数字の集合として \{ -2, -1, 0, 1, 2 \} を採用しているとしましょう。 そうすると、

 2S_j\sum\limits_{k=j+1}^{\infty}Q_kr^{-k} + \left(\sum\limits_{k=j+1}^{\infty}Q_kr^{-k}\right)^2 \le 2S_j\sum\limits_{k=j+1}^{\infty}2\cdot4^{-k} + \left(\sum\limits_{k=j+1}^{\infty}2\cdot4^{-k}\right)^2 = \frac{4}{3}S_j4^{-j}+\frac{4}{9}4^{-2j}  2S_j\sum\limits_{k=j+1}^{\infty}Q_kr^{-k} + \left(\sum\limits_{k=j+1}^{\infty}Q_kr^{-k}\right)^2 \ge 2S_j\sum\limits_{k=j+1}^{\infty}(-2)\cdot4^{-k} + \left(\sum\limits_{k=j+1}^{\infty}-2\cdot4^{-k}\right)^2 = -\frac{4}{3}S_j4^{-j}+\frac{4}{9}4^{-2j}

という上限と下限が求まります。

ここからPDプロットの境界線の表式を求めます。  Q_{j+1}は、以下の式を満たすように選択しなければいけません。

 -\frac{4}{3}S_{j+1}4^{-j-1}+\frac{4}{9}4^{-2j-2} \le w_{j+1} \le \frac{4}{3}S_{j+1}4^{-j-1}+\frac{4}{9}4^{-2j-2}

ここで、 S_{j+1} = S_j + Q_{-j-1}r^{j+1} w_{j+1} = w_j - 2S_jQ_{j+1}4^{-j-1}-(Q_{j+1})^2{}4^{-2j-2} を使って分解して、 Q_{j+1}が式に陽に含まれるようにします。式の左側と右側はほぼ同じなので左側だけ計算しましょう。

 -\frac{4}{3}S_j4^{-j-1}-\frac{4}{3}Q_{j+1}4^{-2j-2}+\frac{4}{9}4^{-2j-2} \le w_j - 2S_jQ_{j+1}4^{-j-1}-(Q_{j+1})^2{}4^{-2j-2}

移項と w_j = R_{j+1}r^{-j-1}の代入を行うと、以下のようになります。

 -\frac{4}{3}S_j4^{-j-1}-\frac{4}{3}Q_{j+1}4^{-2j-2}+\frac{4}{9}4^{-2j-2}+2S_jQ_{j+1}4^{-j-1}+(Q_{j+1})^2{}4^{-2j-2} \le R_{j+1}4^{-j-1}

全体を 4^{j+1}倍します。

 -\frac{4}{3}S_j-\frac{4}{3}Q_{j+1}4^{-j-1}+\frac{4}{9}4^{-j-1}+2S_jQ_{j+1}+(Q_{j+1})^2{}4^{-j-1} \le R_{j+1}

長除法における Dには、開平法における 2S_jが対応していました。 そのため、平方根用のPDプロットの横軸には S_jがとられます。 そこで S_jでくくると、以下のようになります。

 2S_j\left(Q_{j+1}-\frac{2}{3}\right) + \left( (Q_{j+1})^2-\frac{4}{3}Q_{j+1}+\frac{4}{9}\right)4^{-j-1} \le R_{j+1}

この式は、参考文献[1]に記されている式を本稿における表記に変換した*2次の式と一致します。  2S_j\left(Q_{j+1}-\frac{2}{3}\right) + \left(Q_{j+1}-\frac{2}{3}\right)^2{}4^{-j-1} \le R_{j+1}

PDプロットの境界線の表式が求まったため、SRT法で使うテーブルを設計することができます。 このまま設計してもよいのですが、愚直に作ると無駄にSRTテーブルが複雑になってしまいます。 そこで、天下り的ですが、まず私の設計したSRTテーブルを見せ、その中に潜む技巧的な部分を解説していく、という方式にします。

SRT法で浮動小数点数平方根を求めるコード

以下は、浮動小数点数平方根を最近接丸めで求める完全なコードです*3

#include <cassert>
#include <cstdint>

int srt_table( int32_t rem, int32_t div ) {
    assert( 0 <= div && div <= 8 );
    int32_t a1[] = { 6, 7, 8, 8, 9, 10, 10, 11, 11 };
    int32_t a2[] = { 2, 2, 3, 3, 3,  3,  4,  4,  4 };
    if( rem < -a1[div] ) { return -2; }
    if( rem < -a2[div] ) { return -1; }
    if( rem <  a2[div] ) { return 0; }
    if( rem <  a1[div] ) { return 1; }
    return 2;
}

float my_sqrt( float x ) {
    int32_t X = std::bit_cast<int32_t>(x);
    int32_t sign = X >> 31;
    int32_t expo = (X >> 23) & 0xff;
    int32_t mant = X & 0x7fffff;
    if( expo == 0xff && mant != 0 ) { return std::bit_cast<float>(X | 0x00400000); } // sqrt(NaN) = qNaN
    if( expo == 0x00 && mant == 0 ) { return x; } // sqrt(±0.0) = ±0.0
    if( sign ) { return std::bit_cast<float>(0xffc00000); } // sqrt(-X) = qNaN
    if( expo == 0xff && mant == 0 ) { return x; } // sqrt(+Inf) = +Inf
    if( expo == 0 ) { // Subnormal handling
        while( mant <<= 1, (mant & 0x800000) == 0 ) {
            --expo;
        }
    }
    int32_t mx = expo & 1 ? (0x800000 | mant) << 1 : (0x800000 | mant) << 2; // 2 * x
    int32_t quo = expo & 1 ? 0x1600000 : 0x1800000; // 1.375 or 1.5
    int32_t rem = mx - (expo & 1 ? 0x1e40000 : 0x2400000); // 2 * (x - quo * quo)
    for( int i = 1; i < 13; ++i ) {
        int q = srt_table( rem >> 21, (quo >> 21) - 8 );
        rem = (rem << 2) - (q * quo << 1) - (q * q << (24-i*2));
        quo = (quo + (q << (24-i*2))) & 0x3ffffff;
    }
    // Here, quo has <1/3ULP error.
    if( quo & 1 ) {
        quo = rem < 0 ? quo - 1 : quo + 1; // round nearest (never tie)
    }
    int32_t res_expo = (expo >> 1) + (expo & 1) + 63;
    int32_t res_mant = (quo >> 1) & 0x7fffff;
    return std::bit_cast<float>(res_expo << 23 | res_mant);
}

ここで使われているSRTテーブルは、除算に使うこともできます*4。 SRTテーブル部分だけでなく、remquoを更新する部分や丸め部分も除算と共用できます。 したがって、プロセッサに除算器と平方根演算器を両方搭載するのではなく、「除算&平方根演算器」を一つだけ搭載すれば、回路面積の節約になります。

ところで、テーブルの値には一部、上記で求めた境界線を越えている部分があります。 なぜこれが問題にならないのかを、以下で説明していきます。

SRTテーブルの設計

図1: 平方根演算器用のPDプロット
図1は、 j=1,2,\inftyの境界線を書き込んだPDプロットです。 横軸は S_{j-1}、縦軸は R_jです。  j=\inftyの境界線の上に j=2の境界線があり、さらにその上に j=1の境界線があることが分かります。  j=3以降の境界線は書いていませんが、 j=2の境界線と j=\inftyの境界線の間に来ます。

これらの境界線は完全に重なることはないものの、 jごとに表を完全に作り直さないといけないほどはずれていません。 言い換えると、 jの値によらず、同じ表を使うことができます。 その場合、一番厳しい境界線に合わせてテーブルを構成する必要があります*5

図2: 平方根演算器用のSRTテーブルの構成法
図2は、 j=1,2,\inftyの境界線に加え、上記の設計におけるSRTテーブルの値での色分けを書き込んだPDプロットです*6。 「一番厳しい境界線に合わせてテーブルを構成する必要がある」といいながら、いくつかの個所では色分けが境界線を越えてしまっていることが分かります。 具体的には、図3に示した部分は色分けが境界線を越えてしまっています。 なぜこれは問題にならないのでしょうか。
図3: SRTテーブルの色分けが境界線を越えている部分

境界越えα

三か所のαでは、 j=\inftyの境界線がちょうど格子点を通っています。 したがって、有限の jの場合、 S_jが1.125より少し小さい部分、1.5より少し小さい部分、1.875より少し小さい部分において色分けが境界線を越えてしまっています。 しかしこれは問題になりません。 なぜなら、境界線を越えた部分(たとえば、1.125より少し小さい部分)を使うことはないからです。 例えば j=3の場合、  S_jは1,1.125.1.25,1.375,1.5,1.625,1.75,1.875,2.0しかとらないため、1.125より少し小さい値、1.5より少し小さい値、1.875より少し小さい値、になることはありません。

1.125より少し小さい部分での問題を詳しく見ていきます(1.5より少し小さい部分と1.875より少し小さい部分も全く同じ状況です)。  S_jとして取りうる値のうち1.125より小さく1.125に近いものは、 jが大きくなるにつれ1.125に近づいていきます。  S_jとして取りうる値のうち1.125より小さいものの中での最大は、 \frac{9}{8}-2^{-j}と表されます( j \ge 3の場合)。 例えば、 j=4の時1.0625、 j=5の時1.09375、 j=6の時1.109375、などです。 しかしながら、境界線とP=3.0の交点はそれを上回る速度で1.125に近づくため、やはり境界線を越えた部分を使うことはありません。 境界線とP=3.0の交点は、 \frac{9}{8}-\frac{2}{3}\cdot4^{-j}と表されます。 具体的には、 j=4で1.122396、 j=5で1.124349、 j=6で1.124837です。 このような理由から、三か所のαでは見た目上色分けが境界を越えていますが、そのような部分が使われることはなく、問題になりません。

境界越えβ

βでは、 S_jが1.5の時、色分けが j=1の境界線を越えています(D=1.5との交点は -4-\frac{11}{36}であり、少しだけ-4.5より大きいです)。 しかし、この部分は使われることがないので問題になりません。 この部分が使われないのは、入力 N 1\le N \lt 4の範囲に正規化したとき、 1 \le N \lt 2の場合は S_0として1.375を設定しているからです。 SRTアルゴリズムを動かす場合、初期値 S_0は真の解から 2/3までずれた値が許されます。 それにもかかわらず、初期値を一律1.5とせず、 N \lt 2か否かで場合分けしている理由は、この部分を避けることにあります。

境界越えγ

γでは、 S_jが1.25の時、色分けが j=2の境界線を越えています(D=1.25との交点は -3-\frac{143}{144}であり、わずかに-4.0より大きいです)。 しかし、この部分は使われることがないので問題になりません。 この部分が使われないのは、入力 N 1\le N \lt 4の範囲に正規化したとき、 1 \le N \lt 2の場合は S_0として1.375を設定しているからです。 境界越えβを避けるためであれば1.25などを設定してもいいのですが、そうすると S_1が1.25になってしまうことがあります。 そこで初期値を1.375とすることで、 S_1が1.125か1.375となり1.25にならないようにしています。 なお、ここで1.125を設定してしまうのは、 S_1が0.875になることがあってSRTテーブルの範囲をはみ出すので不適です。

その他の境界越え

その他、色分けが-1の下限( j=1)の境界線を越えている場所は無数にありますが、 S_0としてあり得るのは1.375または1.5のみであるため、問題になりません。

除算と違い、D=2.0に対応するエントリが必要

色分けが境界を越えている話とは完全に別の話題になります。 先ほど「SRTテーブルは除算用に使うことができる」と言いましたが、正確には、D=2.0に対応するエントリが必要という点で、除算に使うSRTテーブルと異なっています(D=2.0に対応するエントリ以外は除算と同じであり、テーブルが共用できるということ自体は正しいです)。 除算を行う場合、除数 D 1\le D \lt 2の範囲に正規化できるので、SRTテーブルにD=2.0のエントリは不要でした。 平方根を求める場合、暫定商 S_jは正規化しても 1 \le S_j \le 2の範囲を取るのでSRTテーブルにD=2.0のエントリが必要になります。 一応、初期商として 1.11111..._{(4)} \sim 4/3を採用すれば S_j \lt 2とできるのですが、部分剰余の下位ビットが0ではなくなってしまって面倒です。 そもそも、論理圧縮した回路で実装する場合は、D=2.0に対応するエントリが必要なことはあまり問題にならないものと思われます(テーブルをROMで作る場合は入力が1bit増えてしまって嫌な感じですが……)。

参考文献[2]との比較

参考文献[1]の図25に平方根用のPDプロットが掲載されていますが、残念ながらリンク切れで見ることができません。 「連載 コンピュータアーキテクチャの話」を書籍化した、参考文献[2]を確認してみると、以下の図が掲載されていました。

参考文献[2]のp.201より引用。

……もうめちゃくちゃです。

第一に、凡例があっていません。話の流れからすると、実線は U_1 Min、一点鎖線は L_2 Max、点線は  U_0 Min、破線は L_1 Maxであるべきです。

第二に、線がすべて j=\inftyの物になっています。これは、線とs[j]の交点における縦軸の値がいずれもぴったり0.5の倍数になっていることからわかります。 j=1の線はもっと上に来ます。

第三に、この分割ではSRTテーブルを構成することができません。(そもそも境界線が間違っているのですが、この図の境界線を信じたとしても、)一つのマスに部分商が1でなければいけない部分と2でなければいけない部分が混在しています(具体的には、  0.5 \le s[j\rbrack \lt 0.625, 1.5 \le 4W[j\rbrack \lt 1.75のマスです)。それにも関わらず、この図の下の説明では「この図では,s(j)を0.125刻み,4*W(j)を0.25刻みで目盛り線を書いているが,1つのマスに"0"と"1",あるいは"1"と"2"のように異なる値が入っていることはなく,"0 or 1"あるいは,"1 or 2"とくっつければ,マスごとに選択する s_{j+1}の値を一意に決めることができる。」と書かれています。これは明らかに間違いです。

第四に、上半分しか示されていません。除算器用のSRTテーブルを構成する場合、その境界線は上半分と下半分は対称なので上半分だけ見ていてもよいですが、平方根演算器用のSRTテーブルを構成する場合、その境界線は上半分と下半分で非対称である(有限の jで境界線がずれる方向は、どちらでも上方向)ため、上半分だけを示すのでは不十分です。

めちゃくちゃですが、SRT法で平方根が計算できることを知れたのは参考文献[1]のおかげなので、むしろ感謝したいと思います。

参考文献

*1:手で行う開平法では、与数を2桁ずつ降ろして部分剰余を作り、その部分剰余をもとに部分商を決定しますが、基数は100ではありません。2桁ずつ降ろすのは、暫定商 S_i(上 i+1桁しか数字が入っていなく、その下は0で埋められている)や“余計な項” (Q_i)^2r^{-i}の位置合わせを意識せずに取り扱うための工夫です。最初から与数全てを部分剰余としても、0-0になる桁が多数発生するだけで、アルゴリズムは全く同じように動作します。

*2:本稿では部分商 Q_iを決定する際に参照される部分剰余を R_iとしていますが、参考文献[1]では部分商 s_jを取った後に算出される部分剰余を W(j)としています。そのため、部分剰余の添え字は本稿の方が参考文献[1]よりも1だけ大きいです。

*3:入力としてあり得る全232通りでx86マシンと結果が一致することを確認済みです。一部、処理系規定動作を前提にしています(左シフトが算術シフトであることを期待しています)。

*4:前回示したSRTテーブルとは一部異なりますが、任意性の範囲内です。また、縦軸のスケールが二倍だけ違います(除数が Dであるところが 2S_{j-1}になったことに対応しています)が、これはremをシフトすることで解決可能です。

*5:正確には、何の仮定もなしにはSRTテーブルを構成することはできません。なぜなら、 1\le D \lt 25/24において「-2の上限( j=\infty)」と「-1の下限( j=1)」の上下が入れ替わっており、この部分で選ぶべき部分商は jの値によらずには決定できないからです。単一のテーブルを使うためには、 1\le D \lt 25/24かつ j=1が発生しないような初期商選択を行うなど、SRTループ外の工夫が必要になります。本記事で紹介しているコードも、そのような工夫が行われています。

*6:この図は、PowerPointのグリッド機能を使って作った色分け画像を、Excelのグラフの背景画像として読み込むことで作りました。

SRT除算について

SRT除算は、回復法や非回復法と同様、商を上の桁から求めていくタイプの除算アルゴリズムです。 この時、回復法や非回復法と違い、商を冗長表現で求める点が特徴です。 商を冗長表現で求めることは、各桁の商の選択に多少の自由度を生み、これが回路的なメリットにつながります。

冗長表現とは、一つの数を表す方法が複数あるような記数法のことです。 N進法では通常、各桁の数字はN種類です。 しかしここでもっと多くの種類の数字を使えるようにすると、同じ数を表す方法が複数生まれます。 例えば、二進法で各桁の数字を0, 1, 2, 3から選ぶとすると、十三は1101とも1013とも213とも表せます(ほかにもいろいろあります)。

図1: 長除法(割り算の筆算)における部分商と部分剰余

商を冗長表現で求めることのメリットを考える前に、まずは普通の、冗長でない表現で商を求める場合を考えます。 冗長でない表現で商を求めるアルゴリズムである、回復法や非回復法を例にとってみます。 回復法や非回復法では、その時点での余り(部分剰余)を用いて各桁の商(部分商)を上の桁から順に決定していきます。 この時、部分商の決定を間違えることは許されません。 したがって、以下のように部分商を決定する際に部分剰余の正確な値が必要となります*1

  • 回復法では、除数と部分剰余を比べ、除数≧部分剰余なら部分商として1を立て、除数<部分剰余なら部分商として0を立てる。
  • 非回復法では、部分剰余の符号ビットが0なら部分商として1を立て、部分剰余の符号ビットが1なら部分商として-1を立てる。

商を冗長表現で求める場合、部分商を決定する際のある種の“間違い”を許容できます。 例えば、最終的な答えが十三の場合を考えてみます。 冗長でない二進法で求める場合、部分商の決定に間違いは許されず、順に1, 1, 0, 1と選んで1101とする必要があります。 一方、各桁に0, 1, 2, 3を使える冗長な二進法で求める場合、1, 0, と選んでしまっても許容範囲内です。 なぜなら、その後に2, 1と選んで1021としたり、あるいは1, 3と選んで1013としても、十三になるからです。 このように、特定の桁で多少“間違った”部分商を立ててしまっても、後続の部分商の選択をうまくやればリカバリーが可能です。

部分商を決定する際のある種の“間違い”が許容できることは、以下のような回路的メリットにつながります。

  • 部分商を決定する際、除数と部分剰余を全桁見る必要はなく、除数の上位数ビットと部分剰余の上位数ビットのみを見ればよい
    • テーブル引き*2での実装が可能になる。部分商を二~三桁同時に求めることも可能(Radix-4、Radix-8)*3
    • テーブル引きでなく回復法のように引き算した結果の符号を使って部分商を選択する実装の場合でも、その引き算を誤差あり(上位ビットだけの計算)で行うことができるようになる(参考:[2])
  • 部分剰余を正確に求める必要はなく、その上位数ビットのみが±1程度の誤差で求まればよい
    • 部分剰余は桁上げ保存加算器(と上位数ビットを確定させる小さな桁上げ先見加算器)で求めればよい。キャリーの伝搬を待つ必要がなくなる

さて、どの程度“間違った”部分商を立ててしまってもリカバリーが可能なのでしょうか。 これは、後続の桁の部分商をなるべく大きくした場合となるべく小さくした場合を考えることで求めることができます。

(なるべく大きくした場合)先の記数法では、後続の部分商を33333...とするのが最も大きくなります。 自分より下の桁が33333...となる時、その重みは自分の桁に対して0.33333...になります。 この値は三に等しいので、「真の商より(自分の桁で)三まで小さい商を立ててしまってもリカバリ可能」ということになります。

(なるべく小さくした場合)先の記数法では、後続の部分商を00000...とするのが最も小さくなります。 自分より下の桁が00000...となる時、その重みは自分の桁に対して0.00000...になります。 この値は零に等しいので、「真の商より(自分の桁で)零まで大きい商を立ててしまってもリカバリ可能(=ちょっとでも大きい商を立ててしまうとリカバリ不能)」ということになります。

ただし、これらの猶予は、部分剰余を正確に求めない場合には狭くなるので注意が必要です(PentiumのFDIVバグなど。参考:[2])。

ここまでがSRT除算の基本になります。 以下では、Radix-4・最小冗長表現・テーブル引き・部分剰余を正確に求める場合のSRT除算について詳しく見ていきます。 Radix-4とは四進で(=二進で二桁ずつ)商を求めること、最小冗長表現とはk進法に対して各桁の数字がk+1種類あることを意味します。

Radix-4・最小冗長表現・テーブル引き・部分剰余を正確に求める場合のSRT除算

Radix-4・最小冗長表現のSRT除算では、各桁の数字として-2, -1, 0, 1, 2を使います。 部分剰余を正確に求める場合、リカバリ可能な範囲を考えてみると、

  • 下の桁を22222...とすると最大で、その重みは自分の桁に対して0.22222...=三分の二
  • 下の桁を(-2)(-2)(-2)(-2)(-2)...とすると最小で、その重みは自分の桁に対してマイナス三分の二

となるので、真の商からプラスマイナス三分の二の範囲内に位置する商を立てればリカバリー可能です。 真の商からプラスマイナス二分の一の範囲内で商を立てなければいけなかった非回復法と比較してみてください。

以下では、除数・部分剰余・商の関係を見ていきます。 先ほどのリカバリ可能な範囲の議論から、除数Dと部分剰余R*4について、以下の関係が成り立ちます。

  • (4/3)D ≦ Rである時に限って、商として2を立ててよい。
  • (1/3)D ≦ R ≦ (5/3)Dである時に限って、商として1を立ててよい。
  • (-2/3)D ≦ R ≦ (2/3)Dである時に限って、商として0を立ててよい。
  • (-5/3)D ≦ R ≦ (-1/3)Dである時に限って、商として-1を立ててよい。
  • R ≦ (-4/3)Dである時に限って、商として-2を立ててよい。

ここで注目したいのは、これらの範囲には重複があることです。 例えば、(4/3)D ≦ R ≦ (5/3)Dの範囲の場合、商として1を立てても2を立ててもよいことになります。 こういった部分では、“間違い”が許容されるわけです。 一方、(2/3)D < R < (4/3)Dといった範囲では、“間違い”は許容されず、商として必ず1を立てなければなりません。

この重複の範囲をうまく利用すると、テーブル引きで商を求めることが可能です。 これを考えるために、まずはPDプロット*5を作ります。 横軸に除数D、縦軸に部分剰余Rを取って、先ほどの「(4/3)D ≦ R」などの境界線を引くと、図2のようになります*6。 ただし、除数Dは1≦D<2の範囲のみを図示してあります。 これは、この範囲内に正規化できることを考慮しています。 また、上下対称であるため、上半分のみを示しています。

図2: 最小冗長のRadix-4 SRTにおけるPDプロット

テーブル引きに使う表を作るには、これを格子状に区切り、各マスを商に対応した一つの色で塗ります。 それをやったのが図3のようになります(網掛け部分は、塗る色にいくらかの任意性があります)。 選択される商が許容される範囲内であることを確認してみてください。 境界線は斜めですが、猶予の部分をうまく活用することで、格子であっても許容された範囲内の商を選択することができています。 また、図3では、横軸を1/8刻み、縦軸を1/2刻みで区切ってあります。 除数Dは1≦D<2、部分剰余Rは-16/3<-(8/3)D≦R≦(8/3)D<16/3の範囲内であるため、これは除数Dの小数部の上位3bitと部分剰余Rの上位6bitのみでテーブル引きできることを意味します。

図3: 最小冗長のRadix-4 SRTで使用するテーブル構成法

これを使った浮動小数点数除算器のサンプルコードは以下になります。

int srt_table( int32_t rem, int32_t div ) {
  int d = ( div >> 20 ) - 8;
  assert( 0 <= d && d < 8 );

  int th12[] = { 6, 7, 8, 9, 9, 10, 11, 12 };
  int th01[] = { 2, 2, 3, 3, 3,  3,  4,  4 };

  rem >>= 21;
  if( rem < -th12[d] ) { return -2; }
  if( rem < -th01[d] ) { return -1; }
  if( rem <  th01[d] ) { return  0; }
  if( rem <  th12[d] ) { return  1; }
  return 2;
}

float srt_divider( float x, float y ) {
  assert( 1.0f <= x && x < 2.0f );
  assert( 1.0f <= y && y < 2.0f );

  int32_t X = x * 0x1.p+23;
  int32_t Y = y * 0x1.p+23;

  if( x < y ) { X *= 2; }

  int32_t quo = 0;
  int32_t rem = X;

  for( int i = 0; i < 13; ++i ) {
    const int q = srt_table( rem, Y );
    rem = ( rem - q * Y ) << 2;
    quo = ( quo << 2 ) + q;
  }
  // Here, quo has <1/3ULP error.
  
  // round nearest (never tie)
  if( quo & 1 ) {
    quo = rem < 0 ? quo - 1 : quo + 1;
  }

  float result = quo * ( x < y ? 0x1.p-25 : 0x1.p-24 );
  assert( result == x / y );
  return result;
}

参考文献

*1:非回復法では部分剰余の符号ビットだけ分かれば十分ですが、符号ビットを求めるのには全桁を求める場合とほとんど変わらない回路遅延が発生するので、「部分剰余の正確な値が必要」としました。

*2:テーブル引きと言っても、実用的にはROMを作るのではなく論理圧縮した回路とするようです。

*3:四桁(Intel PenrynのRadix-16 Dividerとか)や六桁(ARMのRadix-64 Floating-Point Dividerとか)を一サイクルで求めるものもありますが、これらは一サイクル内でRadix-4を複数回繰り返すという手法です。

*4:Partial reminderでPの記号を使う文献[1]もありますが、本稿ではより性質を表しているRを選択しました。文献[2]もRを使用しています。

*5:Partial reminderとdivisorでPDプロットというのが一般的[1][2]なようで、本稿でもそれに倣いました。文献[2]は部分剰余にRを使っているにもかかわらずPDプロットと呼んでいるので、RDプロットとは言わないのが普通のようです。

*6:この図は、PowerPointの6グリッドモードと8グリッドモードを駆使して作りました。

CORDICを使ってsin,cos,sinh,coshなどを求める

CORDICアルゴリズムという、乗算なしにsin/cos/sinh/coshを求めるアルゴリズムがあります。 乗算器すらないような非常に小規模なプロセッサでsinやcosが必要になった場合などに有用です。

sin/cosを求める場合

戦略

加法定理を思い出します。

  •  \sin(\alpha + \beta) = \sin\alpha\cos\beta + \cos\alpha\sin\beta
  •  \cos(\alpha + \beta) = \cos\alpha\cos\beta - \sin\alpha\sin\beta

 \cos\betaを因数としてくくりだすと、以下のようになります。

  •  \sin(\alpha + \beta) = \cos\beta(\sin\alpha + \cos\alpha\tan\beta)
  •  \cos(\alpha + \beta) = \cos\beta(\cos\alpha - \sin\alpha\tan\beta)

したがって、「 \tan\beta \cos\betaが分かっていれば、 \sin\alpha \cos\alphaから \sin(\alpha+\beta) \cos(\alpha+\beta)が求められる」ということになります。

 \tan\betaが二のべき乗になっていれば、 \tan\betaを掛けている部分は単純なシフトで済みます。 したがって、 \tan\betaが二のべき乗になるような \betaをうまく選び、これを足したり引いたりすることで求めたい角度に近づけていくという戦略を取ります。

 \cos\betaについては、全体にかかっているので、最後にまとめて計算します。 実は、この「最後にまとめて計算」するときにかけるべき値は、求めたい角度に近づけていくときの経路によらず事前に計算した値とすることができます。 これは、 \cos\beta \betaの符号によらず同じ値になることを利用しています。 求めたい角度に近づけていく際、「 \betaを足す、または足さない」のように求めたい角度に近づけていくのではなく、「 \betaを足す、または -\betaを足す」のように求めたい角度に近づけていきます。 このようにすることで、「最後にまとめて計算」するときにかけるべき値が、どう近づけていったかによらずに決まります。

準備

CORDICアルゴリズムでは、以下の値を事前に計算して持っておきます。

  •  \beta_k = \arctan(1/2^k) (k=0,1,\dots N)
  •  M = \prod_{k=0}^{N}\cos(1/2^k)

計算

CORDICアルゴリズムでは、sin/cosを求めたい角度 \thetaを以下のように分解します。

 \theta = \pm\beta_0\pm\beta_1\pm\dots\pm\beta_N+\varepsilon

ここで、 \pmについては、誤差 \varepsilonが十分小さくなるよう、うまく符号を選びます。 とはいっても、上位桁から貪欲に選んでいく、という方法が実用上十分です。 選んだ符号を s_kとすると、以下のように書けます。

 \theta = s_0\beta_0+s_1\beta_1+\dots+s_N\beta_N+\varepsilon

すると、以下のように漸化式を順に計算していけば、 \sin \theta \cos \thetaが求まります。

  •  x_0 = 0.0
  •  y_0 = 0.0
  •  x_{k+1} = x_k + s_k y_k2^{-k}
  •  y_{k+1} = y_k - s_k x_k2^{-k}

  •  \sin\theta \simeq Mx_N \cos\theta \simeq My_N

コード

まずは説明のために分かりやすいC++コードを紹介します。

std::pair<double, double> cordic_sincos( double a ) {
        double angle = 0.0;
        double x = 1.0;
        double y = 0.0;
        double M = 1.0;

        for( int i = 0; i < 53; ++i ) {
                double t = std::ldexp(1.0, -i);
                double b = std::atan(t);
                if( angle > a ) {
                        std::tie(x, y) = std::tuple(x + y * t, y - x * t);
                        angle += -b;
                        M *= std::cos(-b);
                } else {
                        std::tie(x, y) = std::tuple(x - y * t, y + x * t);
                        angle += b;
                        M *= std::cos(b);
                }
                M *= std::cos(b);
        }
        return { y * M, x * M };
}

このコードでは毎回Mを計算していますが、実際には分岐の方向によらず同じ値になるため、事前計算可能です。 また、xを初期化する際にMで初期化してしまえば、最後の乗算も不要になります。

それらの最適化を行い、また整数演算だけで計算できるようにしたコードが以下になります。

std::pair<double, double> int64_cordic_sincos( double a ) {
        double angle = 0.0;
        std::int64_t X = 0x4dba76d421af2d34; // Πcos(2^-k)
        std::int64_t Y = 0;
        for( int i = 0; i < 53; ++i ) {
                double t = std::ldexp(1.0, -i);
                double b = std::atan(t);
                if( angle > a ) {
                        std::tie(X, Y) = std::tuple(X + (Y>>i), Y - (X>>i));
                        angle -= b;
                } else {
                        std::tie(X, Y) = std::tuple(X - (Y>>i), Y + (X>>i));
                        angle += b;
                }
        }
        return { Y/0x1.p+63, X/0x1.p+63 };
}

このように、確かに加減算とシフトのみで計算できていることが分かります。 CORDIC法は、固定小数点数で求めたい場合には非常に有力な手法です。

sinh/coshを求める場合

事前に計算する値が \cos \arctanではなく \cosh \mathrm{arctanh}になること、加法定理の符号が一部異なることを除いて、ほぼ同じ……ではありません。

 \beta_k = \arctan(1/2^k) (k=0,1,\dots N)はこれらを組み合わせて適当に足し引きすることで入力 \thetaを精度良く表すことができるのですが、 \beta_k = \mathrm{arctanh}(1/2^k) (k=1,2,\dots N)はそうではありません。 具体的には、例えば0.0付近の数を精度良く表すことができません。  \mathrm{arctanh}(1/2) - \sum_{k=2}^{\infty}\mathrm{arctanh}(1/2^k) > 0.043であり、これより絶対値の小さい値を精度良く表すことができません。

この問題に対しては、 \beta_kの複数使用を許した分解を考えればよいです。 例えば、全部二回使って、 \theta = \pm\beta_1\pm\beta_2\pm\beta_2\pm\beta_3\pm\beta_3\pm\dots\pm\beta_N\pm\beta_N+\varepsilonのようにします。 参考文献[2]によれば、重複利用するのは三回に一回程度でよいようです。 これは、 \mathrm{arctanh}(1/2^n) -  \sum_{k=n+1}^{\infty}\mathrm{arctanh}(1/2^k) \lt \mathrm{arctanh}(1/2^{3n+1})が成り立つためです。 具体的なコードは以下のようになります。

std::pair<double, double> cordic_sinhcosh( double a ) {
        double angle = 0.0;
        double x = 1.0;
        double y = 0.0;
        double M = 1.0;

        for( int i = 1; i < 53; ++i ) {
                double t = std::ldexp(1.0, -i);
                double b = std::atanh(t);
                if( angle > a ) {
                        std::tie(x, y) = std::tuple(x - y * t, y - x * t);
                        angle += -b;
                        M *= std::cosh(-b);
                } else {
                        std::tie(x, y) = std::tuple(x + y * t, y + x * t);
                        angle += b;
                        M *= std::cosh(b);
                }
                if( i != 1 && i % 3 == 1 ) {
                        if( angle > a ) {
                                std::tie(x, y) = std::tuple(x - y * t, y - x * t);
                                angle += -b;
                                M *= std::cosh(-b);
                        } else {
                                std::tie(x, y) = std::tuple(x + y * t, y + x * t);
                                angle += b;
                                M *= std::cosh(b);
                        }
                }
        }
        return { y * M, x * M };
}

CORDICを逆に使う場合

CORDICを逆に使うと、 \arctan \mathrm{arctanh}を求めることができます。 逆に使うとは、以下のようなことを意味します。

  • 普通に使うときは、分解結果 s_0\beta_0+s_1\beta_1+\dots+s_N\beta_Nと入力 \thetaの誤差 \varepsilonが0に近くなるように符号 s_kを定めつつ、 x, yを更新していきました。
    • 入力が \theta、出力が x, yです。
    • 初期点(x, y)から \theta回転させたとき、(x, y)がどこにあるかを計算してことになります。
  • 逆に使うときは、 yが0に近くなるように符号 s_kを定めつつ、 x, yを更新していきます。副産物として、分解結果 \theta = s_0\beta_0+s_1\beta_1+\dots+s_N\beta_Nが得られますが、これが求めるものです。
    • 入力が x, y、出力が \thetaです。
    • 初期点(x, y)から(x', 0)に近づけるように回転したとき、何度回転させたのかを計算していることになります。

なお、sin/cosを求めるときのCORDICを逆に使うと \sqrt{x^2+y^2}を、sinh/coshを求めるときのCORDICを逆に使うと、 \sqrt{x^2-y^2}を、それぞれ求めることができます。 この時は、出力は \thetaではなく xになります(点(x, y)を通る円・双曲線のx切片を求めていることになります)。 つまり、 \thetaは全く参照されないので、平方根を求めるだけであれば \arctan \mathrm{arctanh}のテーブルは不要です。 さらに、sinh/coshを求めるときのCORDICに x = b + 1/4, y = b - 1/4を入力して逆に使うことで、 \sqrt bを求めることもできます。

参考文献

浮動小数点数演算における保護桁について

浮動小数点数同士を足したり引いたりする場合、その数学的な結果は浮動小数点数で表せるとは限りません。 そのため、結果を適当に浮動小数点数に戻す、丸めが行われます。 丸めを行うためには、最終的な結果に使われる桁数より少し多く情報を保持しておく必要があります。 そのような部分のことを、保護桁と言います。

正確に丸める最も簡単な方法は、丸める前の結果を正確に求める、というものです。 しかしそうするには、かなり大きなビット幅の加算器が必要になります。 幸運なことに、これを回避する方法があります。 最終的な結果に使われる桁数よりもずっと離れた桁の計算はしないというものです。 具体的には、最終的な結果に使われる桁数に加えて3bit余分にあれば、IEEE754で規定された丸めを行うことができます。

この記事では、なぜ余分に3bit必要なのかを説明します。 この3bitには名前がついていて、それぞれ上の桁から「ガードビット」「ラウンドビット」「スティッキービット」となっています。 ガードビットは保護桁の最上位のビット、ラウンドビットはその次のビットを意味しています。 スティッキービットは特殊で、その次以降が全て0ではない(1がいくつか立っている)ことを意味するビットです。

以下では、説明が冗長になることを防ぐため、以下の太字用語を定義して使います。

結論から書くと、以下のようになります。

  • 加算の結果を偶数丸めする場合、ラウンドビットは不要で、ガードビットとスティッキービットのみが必要
  • 減算の結果を偶数丸めする場合、ガードビットとスティッキービットだけでなくラウンドビットも必要
    • 加算と違ってラウンドビットが必要なのは、繰り下がりが発生してガードビット部分が最終結果に使われる場合があるため

加算する場合

浮動小数点数では符号+指数部+絶対値表現なので、以下では符号のことは考えずに絶対値部分のみを考えます。

加算する場合、以下のように桁合わせをしてから加算器に入力することになります。 ここで、桁合わせしてはみ出た部分は、最終結果には直接的には使われません。 したがって、この部分が保護桁であるということになります。 なお、最上位桁に繰上りが発生した場合は、最終結果に使われない部分が1bit増えますが、保護桁が不足するのではなく余るだけなので、問題にはなりません。

図1: 加算する場合の様子

桁合わせしてはみ出た部分については、片側のオペランドが0なので、特に計算せずガードビット(G)とスティッキービット(S)が求まります。

ガードビットが必要なのは、最近接丸めを実現するためです。 0.5より少ないか多いかを判断するために必要になります。 スティッキービットが必要なのは、偶数丸めを実現するためです。 0.5ちょうどなのか0.5より多いのかを判断するために必要になります。 つまり、四捨五入丸め(round nearest, ties to away)をサポートすればよい場合、スティッキービットは不要です。 また、ゼロへの丸めをサポートすればよい場合、ガードビットもスティッキービットも不要です。 他にも丸めモードは多数ありますが、ガードビットとスティッキービットさえあればいずれにも対応できます。

丸めモードごとの各ビットの要不要は以下のようになります。 ここで、「ガードビット不要・スティッキービット必要」とは、はみ出した部分に1があるかを判定するスティッキービットのみが必要であるという意味です。

丸めモード G S
ゼロへの丸め 不要 不要
正の無限大への丸め 不要 必要
負の無限大への丸め 不要 必要
偶数丸め 必要 必要
奇数丸め 不要 必要
四捨五入丸め 必要 不要

偶数丸めを行う場合、ガードビットとスティッキービットの使い方は以下のようになります。

G S 動作
0 0 そのまま
0 1 そのまま
1 0 結果の仮数部の最下位が1ならば1を加える
1 1 結果の仮数部に1を加える

減算する場合

減算する場合も、以下のように桁合わせをしてから加算器に入力することになります。 ここで、桁合わせしてはみ出た部分が保護桁なのですが、注意するべき点があります。 それは、最上位桁で繰り下がりが発生した場合、最終結果に使われないとは言い切れない点です。 ここでは、以下の二通りに分類して考えてみます。

指数部が近く、多数桁の桁落ちが発生しうる場合

オペランドの指数部の差が0か1である場合が相当します。 この場合、桁落ちが発生しますが、結果は浮動小数点数で正確に表せます。 したがって、桁合わせしてはみ出た部分(高々1bit)を含めて加算器に入力し、出力を桁合わせすれば計算完了です。

指数部が近くないので、多数桁の桁落ちが発生しない場合

オペランドの指数部の差が2以上である場合が相当します。 多数桁の桁落ちが発生しないとは言っても、最上位桁の繰り下がりが発生する場合はあります。 その場合、最終結果の有効数字が1bit足りなくなるので、保護桁から捻出する必要が発生します。 よって、加算の時と比べてもう1bit必要で、ガードビット(G)・ラウンドビット(R)・スティッキービット(R)の3bitが必要です。

図2: 減算する場合の様子

足し算の時はあふれたビットがそのまま保護桁となりましたが、引き算の場合はあふれたビットも含めて引き算をして保護桁の値を決定します。 スティッキービットを決めるためには完全な引き算は必要なく、関係するビットを全てORしたものを減算器に入力すれば十分です。 つまり、3bit分だけ幅広だけの減算器が必要ということになります。

減算の場合、四捨五入すればいい場合でもスティッキービットまで必要です(五捨五超入なら不要です)。

偶数丸めを行う場合、保護桁の使い方は以下のようになります。

最上位桁の繰り下がり G R S 動作
なし 0 x x そのまま
なし 1 0 0 結果の仮数部の最下位が1ならば1を加える
なし 1 0 1 結果の仮数部に1を加える
なし 1 1 x 結果の仮数部に1を加える
あり - 0 x そのまま(Gが仮数部の最下位ビットとなる)
あり - 1 0 Gが1ならば結果の仮数部に1を加える
あり - 1 1 結果の仮数部に1を加える

なお、加算の時と回路を共有することを考えると、桁合わせの際に1bitずらして配置することで以下を共通化できるので、そういった構成法がよいかもしれません。

  • 桁合わせに必要な1bit右シフト回路
    • 1bitずらさないと、1bit左シフト回路も必要になる
  • ラウンドビットが不要
  • 丸める際の補正が0か+1しか出てこない

参考文献

qsortの比較関数の書き方について

C言語標準ライブラリに含まれるqsortは、配列のソートを行ってくれる関数です。 第一引数にソート対象の配列、第二引数に配列の要素数、第三引数一要素のサイズ、第四引数に比較関数、のように渡すことで、ソートを行ってくれます。 ここで、第四引数に渡す比較関数をどのように書くかについてまとめてみます。

比較関数の仕様

int comp( const void* lhs, const void* rhs )のような関数を書く必要があります。ここで、返り値は、

  • lhsの方が小さい(より正確には、「配列の先頭に近い方に配置してほしい」)場合は、負の値
  • lhsrhsが等しい場合は、0
  • lhsの方が大きい(より正確には、「配列の末尾に近い方に配置してほしい」)場合は、正の値

を返す必要があります。

C++におけるstd::sortの比較関数のように二値を返すのではなく、三値(以上)を返す必要がある点に注意します。

比較対象がcharの場合

比較対象がcharの場合、うまい方法があります。 例えば昇順にソートしたい場合、以下のように書けば、仕様を満たせます。

int comp( const void* lhs, const void* rhs ) {
    return *(char*)lhs - *(char*)rhs;
}

ここで、*(char*)lhs*(char*)rhschar型ですが、算術演算のオペランドになっているため、int型に変換されてから演算が行われます。 したがって、この計算でオーバーフローは発生せず、仕様通り「lhsの方が小さい時は負の値、等しければ0、lhsの方が大きい場合は正の値」を返すことになります。

比較対象がintの場合

比較対象がintの場合、先のような引き算を使った方法はうまくいきません。 なぜなら、int同士の引き算の結果は、intに収まる保証がないためです。 long longなど、int型より大きい型を使えばよいと思うかもしれませんが、返り値をintに収める必要があるため、やはりうまくいきません。

※あまりにも絶対値の大きな値が来ないことが分かっている場合など、オーバーフローが発生しないことが確実な場合、引き算を使った方法を使っても安全です。

このような場合、以下のような分岐を用いた比較関数を書くことが紹介されることがあります。

int comp( const void* lhs, const void* rhs ) {
    int L = *(int*)lhs;
    int R = *(int*)rhs;
    if( L < R ) { return -1; }
    if( L > R ) { return 1; }
    return 0;
}

しかし、これよりも以下のように書いたほうが高速です。

int comp( const void* lhs, const void* rhs ) {
    int L = *(int*)lhs;
    int R = *(int*)rhs;
    int lt = L < R;
    int gt = L > R;
    return gt - lt;
}

この書き方が高速な理由は、分岐がなくなっているからです。 ソート関数中の分岐を先取りしているだけなので分岐予測ミスの回数自体は変わらなさそうに思いますが、実際には後者の方が分岐予測ミス数が少なくなります。 これはソート関数中の分岐は非常に予測困難であることが関係しているのかもしれません。

測定

分岐を使って比較結果を返すのではなく、return (L > R) - (L < R);と書いたほうが速いのかを実際に確かめてみます。 測定環境は以下です。

  • Ubuntu 20.04.4 LTS (WSL2 on Windows 11)
  • gcc 9.4.0、最適化オプションは-O3 -march=native
  • Intel Core i9 12900K (4880MHzくらいで動いている)

実験に用いたソースコードは以下です。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int comp( const void* lhs, const void* rhs ) {
  long L = *(long*)lhs;
  long R = *(long*)rhs;
  int gt = L > R;
  int lt = L < R;
  return gt - lt;
}

#define SIZE 100000000
long arr[SIZE];

int main() {
        for( int i = 0; i < SIZE; ++i ) {
                arr[i] = rand() + rand() * 0x100000000;
        }
        clock_t start = clock();

        qsort( arr, SIZE, sizeof(long), &comp );

        clock_t end = clock();

        printf( "Elapsed: %f\n", (double)(end - start) / CLOCKS_PER_SEC );
}
int comp( const void* lhs, const void* rhs ) {
  long L = *(long*)lhs;
  long R = *(long*)rhs;
  if (L > R) return 1;
  if (L < R) return -1;
  return 0;
}

測定の結果、前者は10.9±0.1秒程度でソートが終わりましたが、後者は11.3±0.1秒程度の時間がかかりました。 このように、qsortに渡す比較関数中での分岐を減らすことは高速化につながることが分かります。

まとめ

  • charを比較するときには引き算を使うテクニックがある
  • int以上を比較したい場合には引き算のテクニックは使えない(オーバーフローするため)
  • オーバーフローを考慮に入れないといけない場合、(L > R) - (L < R)のように記述して分岐を減らす方が高速