SRT法で平方根を計算する

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

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

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

部分剰余 $R_i$ と除数 $D$ から、部分商 $Q_i$ を何らかの方法で求める。
部分剰余をで更新する
- 被除数 $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$ を何らかの方法で求める。
部分剰余をで更新する。
- 与数 $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は、 $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は、 $j=1,2,\infty$ の境界線に加え、上記の設計におけるSRTテーブルの値での色分けを書き込んだPDプロットです*6。「一番厳しい境界線に合わせてテーブルを構成する必要がある」といいながら、いくつかの個所では色分けが境界線を越えてしまっていることが分かります。具体的には、図3に示した部分は色分けが境界線を越えてしまっています。なぜこれは問題にならないのでしょうか。

境界越えα

三か所のαでは、 $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]を確認してみると、以下の図が掲載されていました。

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

第一に、凡例があっていません。話の流れからすると、実線は $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] Hisa Ando, コンピュータアーキテクチャの話(101) 浮動小数点除算器と平方根演算器 https://news.mynavi.jp/techplus/article/architecture-101/
[2] Hisa Ando, コンピュータ設計の基礎, 毎日コミュニケーションズ (2010)

*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:入力としてあり得る全2³²通りで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のグラフの背景画像として読み込むことで作りました。