平方根を求める回路 - Sis puella magica!

大学の冬休みの課題(ハードウェア構成法講義 - 平木研究室)で25bitの小数P( $1 \le P < 4$ )の平方根A(25bit)を求める回路を書くという課題が出たのでその方法を考えた。

気の向くうちにと勢いで書いた記事なので分かりにくかったり間違えてたりするかも。

問題の設定

小数の平方根も実際にやることは整数の平方根を求めることと同じなので以下では非負整数の平方根を考える。
以下では非負整数のみ考えるので整数といえば非負整数を指すことにする。

25bitという具体的な数ではかえって本質が分かりにくくなるので整数Nと置き換え、問題の設定を(少し改変して)述べると
「2Nビットの整数Pの上位Nビットが与えられる。下位Nビットは0である。また、Pは
$2^{2N-2} \le P < 2^{2N}$
を満たす。
その平方根A(Nビット)を求めよ。」
となる。

丸め処理については1ビット多めに計算することで対応できるのでここでは考えないことにして、 $A^2 \le P$ を満たす最大の整数Aを求める。

mがnビット整数であることは $0 \le m < 2^n$ と同値であり、mの下位nビットが0であることはある整数m'が存在して $m = m' 2^n$ となることと同値であることを注意しておく。

また、nビット整数Iの二進数表現の $2^k$ の位をI(k)と表記する。これはVHDLのstd_logic_vector(n-1 downto 0)を意識した表記である。例えばLSBはI(0)、MSBはI(n-1)となる。また、 $2^i$ の位から $2^j$ の位( $i>j$ )までを取り出したものをI(i..j)のように表記する。これはVHDLではI(i downto j)に対応する。

1.素朴な方法

Aの取りうる範囲に対して二分探索を行うことで平方根を求めることができる。
つまり、 $2^{N-1} \le A < 2^N$ は既に分かっていることに注意すると擬似コードは

A <- 2^(N-1)
for k in N-2..0
  if (A+2^k)^2 <= P
    A <- A+2^k

とすることでAが求まる
この方法では2乗を直接計算しているのであまり良くない

2.開平法

$(A+2^k)^2 = A^2 + 2^{k+1}A + 2^{2k}$
より、1での2乗計算は加算だけで済ませることができる(2のべき乗はシフト演算に過ぎない)
変数Qに $P-A^2$ を保存することにすると、擬似コードは

A <- 2^(N-1)
Q <- P - 2^(2*(N-1))
for k in N-2..0
  if A*2^(k+1) + 2^(2k) <= Q
    A <- A+2^k
    Q <- Q - (A*2^(k+1) + 2^(2k))

これは2進数での開平法を計算していることになる

3.開平法(改良版)

2の方法のQはそのままでは2Nビットの整数として扱うことになるが実はこの中で1が立っているビットはある範囲にかたまっていて、(ループしている間に範囲がスライドしていくが)N+2ビットで表現できる。また、Qに合わせてifの条件文やQの更新の式を書き換えると、シフト演算も1ビットシフトだけでできるようになる。

基本的には2のコードを基本にするのだが、話がややこしくなるので変数の意味とビット幅をまとめておく。

k : ループのカウンタ。n-2から0まで1ずつ減る。
P : 2Nビット。入力された整数(上位Nビット)と0(下位Nビット)。
A : Nビット。最終的にはPの平方根がここに計算される。ループのブロックの先頭では下位k+1ビット(=A(k..0))は0でループの中でA(k)が決定される。初期値は $2^{N-1}$
Q : 2Nビット。 $P - A^2$ に等しい。

このときループのブロックの先頭ではQはQ(N+k+1..k+1)の範囲だけで表現できる(つまり、ループのブロックの末尾ではQ(N+k..k)の範囲で表現される)。qをQ(N+k+1..k+1)とすると

q : N+1ビット。ループのブロックの先頭で $Q = q 2^{k+1}$ が成り立つ。

2のコードをqを用いて書き換えることで、コードが簡略化される。

また、ループのブロックの末尾では、上の条件でkを1減らしたものをそれぞれの変数が満たしている必要がある。つまり、ループ末尾で

A : 下位kビット(=A(k-1..0))は0。
Q : $P - A^2$ に等しい。
q : $Q = q 2^{k}$ を満たす。

A、Qがこの条件を満たすことは明らかなので以下ではqについて考えていく。

QがQ(N+k+1..k+1)の範囲だけで表現できること、つまりqの存在については

ループ先頭で $0 \le Q < 2^{N+N-2+2}$ かつQの下位N-2+1ビットが0
ループ末尾で $0 \le Q < 2^{N+k+1}$ かつQの下位kビットが0

を示せば良い。

最初の方は初期値についての条件だが、これを満たすことは容易に確認できる。

2つ目の方はループ末尾でのA, Qの満たす条件を用いて示す。
Aの下位kビットは0なのでA'(N-kビット整数)が存在して $A = A' 2^k$ を満たす。
$Q = P - A^2 = P - A'^2 2^{2k}$ であり、Pの下位Nビットは0、 $A'^2 2^{2k}$ の下位2kビットは0なのでQの下位 $\min \{N, 2k \}$ ビットは0である。 $0 \le k < N-2$ より $\min \{N, 2k \} \ge k$ なのでQの下位kビットは0である。
また、Aの計算方法を考えると $(A' 2^k)^2 \le P < ((A'+1) 2^k)^2$ が成り立っているので
$0 \le Q = P - (A' 2^k)^2 < (2A'+1) 2^{2k} < 2^{N-k+1} 2^{2k} = 2^{N+k+1}$
よって、 $0 \le Q < 2^{N+k+1}$ かつQの下位kビットが0が示された。

それでは2のコードをqを用いて書き換えていこう。
qの初期値はQ(2N-1..N-1)、つまり $P - 2^{2(N-1)}$ の上位N+1ビット分を取れば良い
ifの条件式 $A 2^{k+1} + 2^{2k} \le Q$ は $Q = q 2 ^ {k+1}$ より $2A + 2^k \le 2q$ と書き換える。両辺ともN+2ビットの整数として計算できる。
Qの更新は
$q 2^k \leftarrow q 2^{k+1} - (A 2^{k+1} + 2^{2k})$
と書き換えて(左辺のQは $q 2^k$ に、右辺のQは $q 2^{k+1}$ になっていることに注意)
$q \leftarrow 2q - (2A + 2^k)$
となる。ここで左辺はN+1ビットだが右辺はN+2ビットで計算したものから最上位1ビット(必ず0になる)を除いたものであることに注意する。
else文を追加しないとif文が不成立の時qが更新されないのでelseに
$q \leftarrow 2q$
も追加する。
改めて擬似コードを書くと

A <- 2^(N-1)
q <- P - 2^(2*(N-1)) の上位N+1ビット
for k in N-2..0
  if 2*A + 2^k <= 2*q
    A <- A+2^k
    q <- 2*q - (2*A + 2^k) の下位N+1ビット
  else
    q <- 2*q の下位N+1ビット

さらに、ループの中に現れるkは2^kの形しかないので、 $c = 2^k$ としてkを除くこともできる。最終的には

A <- 2^(N-1)
q <- P - 2^(2*(N-1)) の上位N+1ビット
c <- 2^(N-2)
while c > 0
  if 2*A + c <= 2*q
    A <- A + c (+ではなくorでも同じ)
    q <- 2*q - (2*A + c) の下位N+1ビット
  else
    q <- 2*q の下位N+1ビット
  c <- (c >> 1) (右に1ビットシフト)

となる。
VHDLのコンパイラがどれほど賢いのか自分はわからないので、最後のkを取り除く変形に意味があるのかは怪しいが、1ビットシフトの方が単純そうな気がする。