除算 (デジタル)

数値的（ディジタル）な除算アルゴリズムはいくつか存在する。それらのアルゴリズムは、低速な除算と高速な除算の2つに分類できる。低速な除算は反復する毎に最終的な商を1桁ずつ生成していくアルゴリズムである。回復型、不実行回復型、非回復型、SRT除算などがある。高速な除算は最初に商の近似値から出発して徐々に正確な値に近づけていくもので、低速な除算よりも反復回数が少なくて済む。ニュートン-ラプソン法とゴールドシュミット法がこれに分類される。

以下の解説では、除算を $Q=N/D$ で表し、

Q = 商 (quotient)
N = 被除数（分子 = numerator）
D = 除数（分母 = denominator）

とする。

余りのある整数除算（符号なし）[編集]

ここで示すアルゴリズムでは、N を D で割って、商 Q と余り R (remainder) を得る。いずれの値も符号なし整数として扱う。

procedure divide_unsigned(N, D: unsigned_integer; var Q, R: unsigned_integer);
const
  n = 32;  (* ビット数 *)
var
  i: integer;
begin
  if D = 0 then raise DivisionByZeroException;

  (* 商と余りをゼロで初期化 *)
  Q := 0;
  R := 0;

  for i := n-1 downto 0 do
  begin
    R := 2 * R;                  (* R を1ビット左シフト *)
    R := R + ((N shr i) mod 2);  (* Rの最下位ビットを被除数のiビット目と等しく設定する *)
    if R >= D then
    begin
      R := R - D;
      Q := Q + (1 shl i);
    end;
  end;
end;

これは、後述の回復型と基本的には同じである。

低速な除算技法[編集]

低速な除算技法は全て次の漸化式に基づいている。

P_{j+1}=R\times P_{j}-q_{n-(j+1)}\times D\,\!

ここで

P_j = 部分的剰余 (partial remainder)
R = 基数 (radix)
q_{n − (j + 1)} = 商のビット位置 n-(j+1) の桁の値。ここでビット位置は最下位ビットを 0、最上位ビットを n − 1 で表す。
n = 商の桁（ビット）数
D = 除数

である。

回復型除算[編集]

回復型（または復元型）除算 (restoring division) を固定小数点数に対して行う場合を解説する。ここで以下を前提とする。

0 ≤ N
0 < D < 1.

商の各桁 q は数字の集合 {0,1} のいずれかである。

二進（基数2）の基本的アルゴリズムは次の通り。

procedure divide_restoring(N: integer_n_bits; D: integer_2n_bits; var Q: integer_n_bits);
const
  n = 32;
var
  i : integer;
  P : integer_2n_bits;
begin
  (* 商をゼロで初期化 *)
  Q := 0;

  (* P と D は N や Q の倍の幅が必要 *)
  P := N;
  D := D shl n;

  for i := n - 1 downto 0 do
  begin
    P := 2 * P - D;        (* シフトした値から減算を試みる *)
    if P >= 0 then
      Q := Q + (1 shl i);  (* 結果ビットは 1 *)
    else
      P := P + D;          (* 新たな部分的剰余は（回復した）シフトした値 *)
  end;
end;

なお、q(i) は商のi番目のビットを意味する。このアルゴリズムでは減算の前にシフトした値 2P をセーブしておいて、回復（復元）させるステップが必要だが、これはレジスタ T を例えば T = P << 1 としておいて、減算 2P − D の結果が負だった場合にレジスタ T を P にコピーすればよい。

不実行回復型除算は回復型除算とよく似ているが、2*P[i] の値をセーブしておく点が異なり、TP[i] ≤ 0 の場合でも D を加算して戻してやる必要がない。

非回復型除算[編集]

非回復型（または非復元型）除算 (non-restoring division) は商の桁の数字として {0,1} ではなく {−1,1} を使用する。引きっ放し法ともいう。二進（基数2）の基本的アルゴリズムは次の通り。

procedure divide_non_restoring(N, D: integer_n_bits; var q: array_n_of_pm1);
const
  n = 32;
var
  i: integer;
  P: integer_n_bits;
begin
  P := N;

  for i := 0 to n - 1 do
  begin
    if P >= 0 then
    begin
      q[(n - 1) - i] := 1;
      P := 2 * P - D;
    end
    else
    begin
      q[(n - 1) - i] := -1;
      P := 2 * P + D;
    end;
  end;
end;

このアルゴリズムで得られる商は、各桁が −1 と +1 であり、通常の形式ではない。そこで通常の二進形式への変換が必要である。例えば、次のようになる。

次の結果を {0,1} の通常の二進形式に変換する :	$Q=111{\bar {1}}1{\bar {1}}1{\bar {1}}$
ステップ:
1. 負の項のマスクを作る:	$N=00010101\,$
2. Nの2の補数を作る:	${\bar {N}}=11101011$
3. 正の項のマスクを作る:	$P=11101010\,$
4. $P\,$ と ${\bar {N}}$ の総和:	$Q=11010101\,$

ここで、 ${\bar {N}}={\rm {{neg}({\rm {{not}(P))=P+1}}}}$ が成り立つことに注意すると、以下のように修正できる。

procedure divide_non_restoring(N, D: integer_n_bits; var Q: integer_n_bits);
const
  n = 32;
var
  i: integer;
  P: integer_n_bits;
begin
  Q := 0;
  P := N;

  for i := 0 to n - 1 do
  begin
    if P >= 0 then
    begin
      Q := Q + (1 shl ((n - 1) - i));
      P := 2 * P - D;
    end
    else
      P := 2 * P + D;
  end;

  Q := 2 * Q + 1;

  if P < 0 then Q := Q - 1;  (* 引き過ぎている場合、戻す *)
end;

SRT除算[編集]

SRT除算の名は、考案者のイニシャル (Sweeney, Robertson, Tocher) に因んだもので、多くのマイクロプロセッサの除算器の実装に使われている。SRT除算は非回復型除算に似ているが、被除数と除数に基づいてルックアップテーブルを使い、商の各桁を決定する。基数を大きくすると一度の反復で複数ビットを求められるため、高速化可能である^[1]。

いわゆるPentium FDIV バグは、このルックアップテーブルの間違いが原因だった。テーブルのエントリのうち、5個のエントリについて、参照されないものであるという誤った前提を立てて設計してしまったため、オペランドの値によってそれを参照するような場合には、結果がごくわずかだがおかしくなった^[2]。

高速な除算技法[編集]

ニュートン-ラプソン除算[編集]

ニュートン-ラプソン除算 (Newton-Raphson Division) は、ニュートン法を用いて $D$ の逆数を求め、その値と $N$ の乗算を行うことで商 $Q$ を求める。

ニュートン-ラプソン除算のステップは次の通り。

除数 ( $D$ ) の逆数の近似値を計算する: $X_{0}$
逆数のより正確な近似値を反復的に計算する: $(X_{1},\ldots ,X_{S})$
被除数と除数の逆数の乗算を行うことで商を計算する: $Q=NX_{S}$

$D$ の逆数をニュートン法で求めるには、 $X=1/D$ で値がゼロとなる関数 $f(X)$ を求める必要がある。明らかにそのようになる関数としては $f(X)=DX-1$ があるが、これは $D$ の逆数が既にわかっていないと使えない。さらに $f(X)$ の高次導関数が存在しないため、反復によって逆数の精度を増すこともできない。実際に使用できる関数は $f(X)=1/X-D$ で、この場合のニュートン-ラプソンの反復は次の式で表される。

X_{i+1}=X_{i}-{f(X_{i}) \over f'(X_{i})}=X_{i}-{1/X_{i}-D \over -1/X_{i}^{2}}=X_{i}+X_{i}(1-DX_{i})=X_{i}(2-DX_{i})

この場合、 $X_{i}$ から乗算と減算だけで計算可能であり、積和演算2回でもよい。

誤差を $\epsilon _{i}=DX_{i}-1\,$ と定義すると

X_{i}={1 \over D}(1+\epsilon _{i})\,

\epsilon _{i+1}=-{\epsilon _{i}}^{2}\,

となる。

除数 D が 0.5 ≤ D ≤ 1 となるようビットシフトを施す。同じビットシフトを被除数 N にも施せば、商は変化しない。すると、ニュートン-ラプソン法の初期値として次のような線形近似が使える。

X_{0}=T_{1}+T_{2}D\approx {\frac {1}{D}}\,

区間 $[0.5,1]$ においてこの近似の誤差の絶対値をなるべく小さくするには、次の式を使用する。

X_{0}={48 \over 17}-{32 \over 17}D\,

^[要出典]

この近似を用いると、初期値の誤差は次のようになる。

\vert \epsilon _{0}\vert \leq {1 \over 17}\approx 0.059\,

この技法の収束は正確に二次的なので、 $P\,$ 桁の二進数の値を計算する場合のステップ数は次のようになる。

S=\left\lceil \log _{2}{\frac {P+1}{\log _{2}17}}\right\rceil \,

ゴールドシュミット除算[編集]

ゴールドシュミット除算の名は Robert Elliott Goldschmidt に因んだもので^[3]、除数と被除数の両方に共通の係数 F_i をかけていき、除数 D が 1 に収束するようにする。すると被除数 N は商 Q に収束する。つまり、以下の式で分母が1になるようにもっていく。

Q={\frac {N}{D}}{\frac {F_{1}}{F_{1}}}{\frac {F_{2}}{F_{2}}}{\frac {F_{\ldots }}{F_{\ldots }}}

ゴールドシュミット除算のステップは次の通り。

乗数となる係数 F_i を推定により生成する。
除数と被除数に F_i をかける。
除数が十分 1 に近くなったら、被除数を返す。さもなくばステップ1に戻ってループする。

0 < D < 1 となるよう N/D を調整済みとし、それぞれの F_i は D から次のように求める。

F_{i+1}=2-D_{i}

除数と被除数にその係数をかけると次のようになる。

{\frac {N_{i+1}}{D_{i+1}}}={\frac {N_{i}}{D_{i}}}{\frac {F_{i+1}}{F_{i+1}}}

k 回の反復で十分なら、 $Q=N_{k}$ となる。

ゴールドシュミット法はAMDの Athlon やその後のモデルで使用されている^[4]^[5]。

二項定理[編集]

ゴールドシュミット法は、二項定理を使ってより単純化した係数を使うことができる。 $D\in ({\tfrac {1}{2}},1]$ となるよう N/D を2の冪でスケーリングすることを前提とする。ここで $D=1-x$ となるよう x を求め、 $F_{i}=1+x^{2^{i}}$ とする。すると次のようになる。

{\frac {N}{1-x}}={\frac {N\cdot (1+x)}{1-x^{2}}}={\frac {N\cdot (1+x)\cdot (1+x^{2})}{1-x^{4}}}={\frac {N\cdot (1+x)\cdot (1+x^{2})\cdot (1+x^{4})}{1-x^{8}}}

$x\in [0,{\tfrac {1}{2}})$ なので、 $n$ ステップ後には $1-x^{2^{n}}$ と 1 の相対誤差は $2^{-n}$ となり、 $2^{n}$ の二進数の精度では 1 と見なせるようになる。このアルゴリズムをIBM方式と呼ぶこともある^[6]。

大きな整数の場合[編集]

ハードウェアの実装に使われている設計技法は、一般に数千桁から数百万桁の十進数値での除算（任意精度演算）には適していない。そのような除算は例えば、RSA暗号の合同式の計算などでよく見られる。大きな整数での効率的除算アルゴリズムは、まず問題をいくつかの乗算に変換し、それに漸近的に効率的な（つまり桁数が大きいほど効率がよい）乗算アルゴリズム（英語版）を適用する。例えば、トム・クック乗算（英語版）やショーンハーゲ・ストラッセン法がある。乗算への変換としては、上述したニュートン法を使った例や^[7]、それより若干高速な Burnikel-Ziegler の除算アルゴリズム（ドイツ語版）^[8] や Barrett reduction アルゴリズムがある^[9]。ニュートン法は同じ除数で複数の被除数に対して除算を行う場合に特に効率的で、除数の逆数を1度計算しておくと、毎回それを流用できる。

定数による除算[編集]

定数を除数とする除算は、その定数の逆数との乗算と等価である。そのため、除数 D がコンパイル時にわかっている場合（定数の場合）、その逆数 (1/D) をコンパイル時に計算すれば、N·(1/D) という乗算のコードを生成すればよいということになる。浮動小数点数の計算であれば、そのまま適用できる。

整数の場合は、一種の固定小数点数による計算に変形する手法がある。まず、算術的に考えてみる。例えば、除数が3の場合、2/3、4/3、256/3などのどれかを使って乗算し、しかる後に2や4や256で除算すればよい。2進法であれば除算はシフトするだけで良い（16ビット×16ビット＝32ビットのような、倍長で演算結果が全て得られる計算機なら、運が良ければ上位16ビットにそのまま解が得られるようにすることもできる）。

これを整数演算でおこなう場合は、256/3は当然正確な整数にはならないので、誤差があらわれる。しかし、シフト幅をより大きくし、値の範囲に注意すれば、常に不正確な部分はビットシフトによって捨てられる^[10]ように変形できることがある。

具体例として32ビットの符号なし整数で、除数が3の場合 $2863311531/2^{33}$ との乗算に変換できる。まず 2863311531 との乗算を行い、その後33ビット右シフトする。この値は正確には $1/2.999999999650754$ である。

場合によっては、定数による除算を一連のシフト操作と加減算に変換できることもある^[11]。特に興味深いのは10による除算で、シフトと加減算で正確な商（と必要なら余り）が得られる^[12]。

脚注[編集]

^ 葛毅、阿部公輝、浜田穂積 (2002-08). “高基数SRT除算の論理回路実現に基づく回路構成と評価”. 情報処理学会論文誌 43 (8). ^{[リンク切れ]}
^ Intel Corporation, 1994, Retrieved 2011-10-19,"Statistical Analysis of Floating Point Flaw"
^ Robert E. Goldschmidt, Applications of Division by Convergence, MSc dissertation, M.I.T., 1964
^ Stuart F. Oberman, "Floating Point Division and Square Root Algorithms and Implementation in the AMD-K7 Microprocessor", in Proc. IEEE Symposium on Computer Arithmetic, pp. 106–115, 1999
^ Peter Soderquist and Miriam Leeser, "Division and Square Root: Choosing the Right Implementation", IEEE Micro, Vol.17 No.4, pp.56–66, July/August 1997
^ Paul Molitor, "Entwurf digitaler Systeme mit VHDL"
^ Hasselström, Karl (2003). Fast Division of Large Integers: A Comparison of Algorithms (PDF) (Master's in Computer Science thesis). Royal Institute of Technology. 2011年3月23日閲覧。
^ Joachim Ziegler, Christoph Burnikel (1998), Fast Recursive Division, Max-Planck-Institut für Informatik, オリジナルの2011-04-26時点におけるアーカイブ。 2021年9月10日閲覧。
^ Paul Barrett (1987). "Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor". Proceedings on Advances in cryptology---CRYPTO '86. London, UK: Springer-Verlag. pp. 311–323. ISBN 0-387-18047-8。
^ Division by Invariant Integers using Multiplication Torbjörn Granlund and Peter L. Montgomery. ACM SIGPLAN Notices Vol 29 Issue 6 (June 1994) 61–72
^ Massmind: "Binary Division by a Constant"
^ R. A. Vowels, "Divide by 10", Australian Computer Journal (24), 1992, pp. 81-85.

外部リンク[編集]

Computer Arithmetic Algorithms JavaScript Simulator – 各種除算アルゴリズムのシミュレータがある。

[1] 葛毅、阿部公輝、浜田穂積 (2002-08). “高基数SRT除算の論理回路実現に基づく回路構成と評価”. 情報処理学会論文誌 43 (8). ^{[リンク切れ]}

[2] Intel Corporation, 1994, Retrieved 2011-10-19,"Statistical Analysis of Floating Point Flaw"

[3] Robert E. Goldschmidt, Applications of Division by Convergence, MSc dissertation, M.I.T., 1964

[4] Stuart F. Oberman, "Floating Point Division and Square Root Algorithms and Implementation in the AMD-K7 Microprocessor", in Proc. IEEE Symposium on Computer Arithmetic, pp. 106–115, 1999

[5] Peter Soderquist and Miriam Leeser, "Division and Square Root: Choosing the Right Implementation", IEEE Micro, Vol.17 No.4, pp.56–66, July/August 1997

[6] Paul Molitor, "Entwurf digitaler Systeme mit VHDL"

[7] Hasselström, Karl (2003). Fast Division of Large Integers: A Comparison of Algorithms (PDF) (Master's in Computer Science thesis). Royal Institute of Technology. 2011年3月23日閲覧。

[8] Joachim Ziegler, Christoph Burnikel (1998), Fast Recursive Division, Max-Planck-Institut für Informatik, オリジナルの2011-04-26時点におけるアーカイブ。 2021年9月10日閲覧。

[9] Paul Barrett (1987). "Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor". Proceedings on Advances in cryptology---CRYPTO '86. London, UK: Springer-Verlag. pp. 311–323. ISBN 0-387-18047-8。

[10] Division by Invariant Integers using Multiplication Torbjörn Granlund and Peter L. Montgomery. ACM SIGPLAN Notices Vol 29 Issue 6 (June 1994) 61–72

[11] Massmind: "Binary Division by a Constant"

[12] R. A. Vowels, "Divide by 10", Australian Computer Journal (24), 1992, pp. 81-85.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]