128bit除算も定数除算最適化したい

一般に整数除算（ここでは剰余演算も含むこととします）を行う命令は遅いです。 そのため、除算はなるべく回避したいものです。

除数がコンパイル時定数の場合、コンパイラの最適化により、除算が取り除かれることがあります。 これを定数除算最適化と呼びます。

定数除算最適化では、除算はいくつかの軽量な命令に分解されます。 二のべき乗で割る場合、単にシフト命令が使われます。 それ以外の数で割る場合、乗算とシフト命令が使われます。 符号付きの除算の場合、零への丸めの部分の取り扱いが面倒ですが、追加のシフト命令と加算命令だけで（分岐命令を使わずに）なんとかできるようです。

さて、この定数除算最適化は、128bit整数の除算の場合は常には行われないようです。 これを何とかしてみたいと思います。 もちろん、手で最適化されたコードを書けばできますが、コンパイル時計算により除数が決まる場合などに対応するのは困難という問題があるなど、柔軟性に欠けます。 そこで、128bit除算を行う、64bit除算のみを用いたC言語ソースコードを作ることにします。 こうすれば、マジックナンバーの導出はコンパイラにまかせることができます。

前提：コンパイラが最適化する場合としない場合

clangは、二のべき乗の時にシフトにする最適化だけ行うようです。 gccも、二のべき乗の時にシフトにする最適化を行います。 それに加えて、除数が小さい時は最適化してくれることが多いようです。 66以下の自然数で割る場合はすべて最適化してくれました。 67,83,101,……で割る場合は__udivti3 （128bit÷128bitを計算する外部ルーチン）の呼び出しになりました。 大きな数でも最適化してくれる数はいくつか存在するようです。 また、二のべき乗以外の時に最適化で出力されるコードは、25～40命令程度です。

実装

基本的な考え方は、64bit÷64bitを基本演算とした長除算を書けばいい、というものです。 これを実装したのが以下のコードになります。

void partial_divide( __uint128_t* R, __uint128_t* Q, unsigned long long D, int shamt ) {
    assert((unsigned long long)(*R >> shamt) == (*R >> shamt));
    *Q += (__uint128_t)((unsigned long long)(*R >> shamt) / D) << shamt;
    *R = *R - (*R >> shamt << shamt) + ((__uint128_t)((unsigned long long)(*R >> shamt) % D) << shamt);
}

void my_divrem( __uint128_t* R, __uint128_t* Q, __uint128_t X, unsigned D ) {
    *R = X;
    *Q = 0;
    partial_divide(R, Q, D, 64);
    partial_divide(R, Q, D, 32);
    partial_divide(R, Q, D, 0);
}

このコードは、Dが32bit整数しか受け付けないという制約がありますが、かなり高速に動作します。

性能測定

実験に使ったコードは以下です。

__uint128_t S = (__uint128_t)1 << 125; // ※グローバル変数にしないと最適化されてしまうので注意
__uint128_t E = S + 1000 * 1000 * 1000;

int main() {
    auto start = std::chrono::system_clock::now();
    __uint128_t sum = 0;
    for( __uint128_t t = S; t < E; ++t ) {
#ifdef G
        sum += t / DD;
#else
        __uint128_t R, Q;
        my_divrem(&R, &Q, t, DD);
        sum += Q;
#endif
    }
    auto end = std::chrono::system_clock::now();
    std::cout << (unsigned long long)sum << " " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << std::endl;
}

私の書いたコードのコンパイルには、clang++12.0.1をおよびg++12.1.0を使いました。コンパイルオプションは-std=c++20 -O3 -march=nativeです。 コンパイルは、AMD EPYC 7702上で行いました（-march=nativeをつけても、以下で使うどのCPU上でも動くコードができました）。

速度の比較相手は、コンパイラの自動最適化が出力した機械語コードです。 これを作るためには、g++12.1.0を使いました。コンパイルオプションは-std=c++20 -O3 -march=native -DGです。

実行にかかった時間は以下のようになりました。

DD=3の場合↓

使ったCPU 自動最適化 my_divrem(clang) my_divrem(gcc)
AMD EPYC 7702 3.0s 3.5s 6.2s
Intel Xeon E5-2643 v3 3.1s 3.4s 5.9s
Intel Core i9-12900KF 1.3s 1.6s 1.9s

DD=67の場合↓

使ったCPU 自動最適化 my_divrem(clang) my_divrem(gcc)
AMD EPYC 7702 27s 3.7s 6.3s
Intel Xeon E5-2643 v3 33s 3.5s 6.3s
Intel Core i9-12900KF 3.9s 1.7s 1.9s

DD=3の場合、my_divrem(clang)は自動最適化よりやや遅い程度となりました。 「自動最適化」はさすがにコンパイラが出力したコードだけあって速いです。 とはいえ、私の強引な実装もそれほど遅くはないようです。 my_divrem(gcc)はかなり遅いですが、これはスピルが大量に発生してしまったのが原因です。 __uint128_tを取り扱うコードはレジスタ圧が高いため、レジスタ割り付けが難しいようです。

DD=67の場合、my_divremはclangもgccも両方とも自動最適化より速くなりました。 「自動最適化」は定数除算最適化をしていないため、除算命令を実行する必要があり低速です。 最近のIntelのCPUでは除算命令が高速化されたのか、除算命令にコンパイルされてもある程度の速さですが、依然として定数除算最適化をした場合よりも二倍以上遅いです。

もっと大きな数での剰余

Dが32bitよりも少し大きいくらいの場合であれば、partial_divideを呼ぶ回数を増やすことで対処可能です。 例えば、Dが42bitに収まる整数であれば、次にようにすればよいです。

void my_divrem( __uint128_t* R, __uint128_t* Q, __uint128_t X, unsigned long long D ) {
    *R = X;
    *Q = 0;
    partial_divide(R, Q, D, 66);
    partial_divide(R, Q, D, 44);
    partial_divide(R, Q, D, 22);
    partial_divide(R, Q, D, 0);
}

つまり、Dがkビットであれば、(64-k)ビットずつずらしながら除算を行えばよいです。 どこまでやればいいかというと、(64-k)*Nが64以上になるところまでです。

剰余環における演算がしたい場合

剰余環における乗算を実装するために128bit整数が必要となっている場合、さらに最適化することができます。 法をMとして、Mが42bitに収まる整数であれば、次のようにして問題ありません。

unsigned long long mul_with_mod( unsigned long long x, unsigned long long y, unsigned long long M ) { // x,yはM未満とする
    __uint128_t Q, R = (__uint128_t)x * y; // Rは84bitに収まる
    partial_divide(&R, &Q, M, 22); // 84-22=64なので64bit除算でOK
    partial_divide(&R, &Q, M, 0); // 上の除算の余りRは42+22=64で64bitに収まっている
    return R;
}

Mがkビットである時、(64-k)ビットずつずらしながら除算をするという点は先ほどと同じです。 一方、どこまでやるかが違っていて、(64-k)*Nが2k-64以上になるところまでです。

例えば、法Mが48bitに収まる整数であれば、partial_divideを三回に増やします。このとき、shamtは32, 16, 0となります。 法Mが51bitに収まる整数であれば、partial_divideを四回に増やします。このとき、shamtは39, 26, 13, 0となります。

まとめ

128bit整数を、32bitに収まる整数定数で割るプログラムを最適化しました。 具体的には、コンパイラが最適化しやすいようなコード（partial_divide）から構成される除算関数（my_divrem）を作りました。 その結果、gccが最適化できないような定数除算（例えば、67で割る）プログラムにおいて除算命令を回避し、大幅に高速化できることが分かりました。 gccが最適化する定数除算（例えば、3で割る）プログラムでも、速度低下はわずかでした。 定数除算最適化を行わないclang向けにコンパイルする場合に利用価値が高そうです。

また、32bitより少し大きいくらい（例えば、42bit）の整数で割る場合にも応用できることを示しました。 さらに、重要な応用として、剰余環における乗算の実装の場合、一般の場合より演算回数を少なくできることを示しました。

cat a.out を防ぐ（その2）

また実行形式ファイルをcatしてしまってコンソールが大変なことになったので、対策を強化しました。

今回はa.outという名前ではなかったので、前回（cat a.out を防ぐ - よーる）の対策が働きませんでした。

どういう文字列がコンソールを破壊するのかを調べてみると、[ESC]cという文字列が原因のようです（こんなに短い文字列で発生すると思っていませんでした。実行ファイルをcatするとほぼ毎回、コンソールが大変なことになるのは、標準ライブラリ由来の文字列が原因だと思っていました）。 参考：ANSI escape code - Wikipedia いろいろ実験してみると、[ESC]とcの間に印字可能文字でないものが挟まっていても同様の動作をするようです。

cat a.outだけに対処するのではなく、根本的にこういった問題に対処するため、.catを以下のように変更しました。

cat $args | sed "s/\x1b//g"

コンソールが壊れる原因は[ESC]（0x1b）なので、それを取り除きます。 こうすることで、実行形式ファイルだけでなく、圧縮ファイルなどに含まれる問題のある文字列にも対処することができるようになりました。

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

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

環境

  • 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テーブル部分だけでなく、remやquoを更新する部分や丸め部分も除算と共用できます。 したがって、プロセッサに除算器と平方根演算器を両方搭載するのではなく、「除算＆平方根演算器」を一つだけ搭載すれば、回路面積の節約になります。

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

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グリッドモードを駆使して作りました。