除算 (デジタル)のソースを表示

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

以下の解説では、除算を <math>Q = N/D</math> で表し、
* ''Q'' = 商 (quotient)
* ''N'' = 被除数（分子 = numerator）
* ''D'' = 除数（分母 = denominator）
とする。

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

<syntaxhighlight lang="pascal">
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;
</syntaxhighlight>

これは、後述の[[#回復型除算|回復型]]と基本的には同じである。

== 低速な除算技法 ==
低速な除算技法は全て次の漸化式に基づいている。
:<math>P_{j+1} = R\times P_j - q_{n-(j+1)}\times D\,\!</math>
ここで
* ''P''<sub>''j''</sub> = 部分的剰余 (partial remainder)
* ''R'' = 基数 (radix)
* ''q''<sub> ''n'' &minus; (''j'' + 1)</sub> = 商のビット位置 ''n-(j+1)'' の桁の値。ここでビット位置は最下位ビットを 0、最上位ビットを ''n''&nbsp;&minus;&nbsp;1 で表す。
* ''n'' = 商の桁（ビット）数
* ''D'' = 除数
である。

=== 回復型除算 ===
回復型（または復元型）除算 (restoring division) を[[固定小数点数]]に対して行う場合を解説する。ここで以下を前提とする。
* 0 &le; ''N''
* 0 < ''D'' < 1.

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

二進（基数2）の基本的アルゴリズムは次の通り。
<syntaxhighlight lang="pascal">
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;
</syntaxhighlight>

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

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

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

<syntaxhighlight lang="pascal">
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;
</syntaxhighlight>

このアルゴリズムで得られる商は、各桁が &minus;1 と +1 であり、通常の形式ではない。そこで通常の二進形式への変換が必要である。例えば、次のようになる。
{| border="0" cellpadding="0"
|次の結果を {0,1} の通常の二進形式に変換する : ||<math>Q = 111\bar{1}1\bar{1}1\bar{1}</math>
|-
|ステップ: ||
|-
|1. 負の項のマスクを作る: ||<math>N = 00010101\,</math>
|-
|2. Nの[[2の補数]]を作る: ||<math>\bar{N} = 11101011</math>
|-
|3. 正の項のマスクを作る: ||<math>P = 11101010\,</math>
|-
|4. <math>P\,</math> と <math>\bar{N}</math> の総和: ||<math>Q = 11010101\,</math>
|}

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

<syntaxhighlight lang="pascal">
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;
</syntaxhighlight>

=== SRT除算 ===
SRT除算の名は、考案者のイニシャル (Sweeney, Robertson, Tocher) に因んだもので、多くの[[マイクロプロセッサ]]の除算器の実装に使われている。SRT除算は非回復型除算に似ているが、被除数と除数に基づいて[[ルックアップテーブル]]を使い、商の各桁を決定する。基数を大きくすると一度の反復で複数ビットを求められるため、高速化可能である<ref>{{Cite journal |url= http://almond.cs.uec.ac.jp/papers/pdf/2002/katsu-ipsj_srt.pdf |title= 高基数SRT除算の論理回路実現に基づく回路構成と評価 |journal=情報処理学会論文誌 |volume=43 |number=8 |date=2002-08| author= 葛毅、阿部公輝、浜田穂積}}{{リンク切れ|date=2020年3月}}</ref>。

いわゆる[[Pentium FDIV バグ]]は、このルックアップテーブルの間違いが原因だった。テーブルのエントリのうち、5個のエントリについて、参照されないものであるという誤った前提を立てて設計してしまったため、オペランドの値によってそれを参照するような場合には、結果がごくわずかだがおかしくなった<ref>Intel Corporation, 1994, Retrieved 2011-10-19,[http://www.intel.com/support/processors/pentium/fdiv/wp/ "Statistical Analysis of Floating Point Flaw"]</ref>。

== 高速な除算技法 ==
=== ニュートン-ラプソン除算 ===
ニュートン-ラプソン除算 (Newton-Raphson Division) は、[[ニュートン法]]を用いて <math>D</math> の[[逆数]]を求め、その値と <math>N</math> の乗算を行うことで商 <math>Q</math> を求める。

ニュートン-ラプソン除算のステップは次の通り。
# 除数 (<math>D</math>) の逆数の近似値を計算する: <math>X_{0}</math>
# 逆数のより正確な近似値を反復的に計算する: <math>(X_{1},\ldots,X_{S})</math>
# 被除数と除数の逆数の乗算を行うことで商を計算する: <math>Q = NX_{S}</math>

<math>D</math> の逆数をニュートン法で求めるには、<math>X=1/D</math> で値がゼロとなる関数 <math>f(X)</math> を求める必要がある。明らかにそのようになる関数としては <math>f(X)=DX-1</math> があるが、これは <math>D</math> の逆数が既にわかっていないと使えない。さらに <math>f(X)</math> の高次[[導関数]]が存在しないため、反復によって逆数の精度を増すこともできない。実際に使用できる関数は <math>f(X)=1/X-D</math> で、この場合のニュートン-ラプソンの反復は次の式で表される。
: <math>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)</math>
この場合、<math>X_i</math> から乗算と減算だけで計算可能であり、[[積和演算]]2回でもよい。

誤差を <math>\epsilon_i = D X_i - 1 \,</math> と定義すると
:<math>X_i = {1 \over D} (1 + \epsilon_i) \,</math>
:<math>\epsilon_{i+1} = - {\epsilon_i}^2 \,</math>
となる。

除数 ''D'' が 0.5&nbsp;≤&nbsp;''D''&nbsp;≤&nbsp;1 となるようビットシフトを施す。同じビットシフトを被除数 ''N'' にも施せば、商は変化しない。すると、ニュートン-ラプソン法の初期値として次のような線形[[近似]]が使える。
:<math>X_0 = T_1 + T_2 D \approx \frac{1}{D} \,</math>

区間 <math>[0.5,1]</math> においてこの近似の誤差の絶対値をなるべく小さくするには、次の式を使用する。
:<math>X_0 = {48 \over 17} - {32 \over 17} D \,</math> {{要出典|date=2012年4月}}

この近似を用いると、初期値の誤差は次のようになる。
:<math>\vert \epsilon_0 \vert \leq {1 \over 17} \approx 0.059 \,</math>

この技法の[[収束]]は正確に二次的なので、<math>P \,</math> 桁の二進数の値を計算する場合のステップ数は次のようになる。
:<math>S = \left \lceil \log_2 \frac{P + 1}{\log_2 17} \right \rceil \,</math>

=== ゴールドシュミット除算 ===
ゴールドシュミット除算の名は Robert Elliott Goldschmidt に因んだもので<ref>Robert E. Goldschmidt, Applications of Division by Convergence, MSc dissertation, M.I.T., 1964</ref>、除数と被除数の両方に共通の係数 ''F''<sub>''i''</sub> をかけていき、除数 ''D'' が 1 に収束するようにする。すると 被除数 ''N'' は商 ''Q'' に収束する。つまり、以下の式で分母が1になるようにもっていく。

:<math>Q = \frac{N}{D} \frac{F_1}{F_1} \frac{F_2}{F_2}  \frac{F_{\ldots}}{F_{\ldots}}</math>

ゴールドシュミット除算のステップは次の通り。
# 乗数となる係数 ''F<sub>i</sub>'' を推定により生成する。
# 除数と被除数に ''F<sub>i</sub>'' をかける。
# 除数が十分 1 に近くなったら、被除数を返す。さもなくばステップ1に戻ってループする。

0&nbsp;<&nbsp;''D''&nbsp;<&nbsp;1 となるよう ''N''/''D'' を調整済みとし、それぞれの ''F<sub>i</sub>'' は ''D'' から次のように求める。

:<math>F_{i+1} = 2 - D_i</math>

除数と被除数にその係数をかけると次のようになる。
:<math>\frac{N_{i+1}}{D_{i+1}} = \frac{N_i}{D_i}\frac{F_{i+1}}{F_{i+1}}</math>

''k'' 回の反復で十分なら、<math>Q=N_k</math> となる。

ゴールドシュミット法は[[アドバンスト・マイクロ・デバイセズ|AMD]]の Athlon やその後のモデルで使用されている<ref>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&ndash;115, 1999</ref><ref>Peter Soderquist and Miriam Leeser, "Division and Square Root: Choosing the Right Implementation", '' IEEE Micro'', Vol.17 No.4, pp.56&ndash;66, July/August 1997</ref>。

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

: <math>
\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}
</math>

<math>x\in[0,\tfrac{1}{2})</math> なので、<math>n</math> ステップ後には <math>1-x^{2^n}</math> と 1 の相対誤差は <math>2^{-n}</math> となり、<math>2^n</math> の二進数の精度では 1 と見なせるようになる。このアルゴリズムをIBM方式と呼ぶこともある<ref>Paul Molitor, [https://mephisto.informatik.uni-halle.de/aufzeichnungen/vhdl/folien/vhdl06_sw.pdf "Entwurf digitaler Systeme mit VHDL"]</ref>。

== 大きな整数の場合 ==
ハードウェアの実装に使われている設計技法は、一般に数千桁から数百万桁の十進数値での除算（[[任意精度演算]]）には適していない。そのような除算は例えば、[[RSA暗号]]の[[合同式]]の計算などでよく見られる。大きな整数での効率的除算アルゴリズムは、まず問題をいくつかの乗算に変換し、それに漸近的に効率的な（つまり桁数が大きいほど効率がよい）{{仮リンク|乗算アルゴリズム|en|multiplication algorithm}}を適用する。例えば、{{仮リンク|トム・クック乗算|en|Toom–Cook multiplication}}や[[ショーンハーゲ・ストラッセン法]]がある。乗算への変換としては、上述した[[ニュートン法]]を使った例や<ref>{{Cite thesis |degree=Master's in Computer Science |title=Fast Division of Large Integers: A Comparison of Algorithms |url= http://www.treskal.com/kalle/exjobb/original-report.pdf |last=Hasselström |first=Karl |year=2003 |publisher=Royal Institute of Technology |accessdate=2011-03-23}}</ref>、それより若干高速な  {{仮リンク|Burnikel-Ziegler の除算アルゴリズム|de|Burnikel-Ziegler-Verfahren}}<ref>{{citation |url=https://domino.mpi-inf.mpg.de/internet/reports.nsf/efc044f1568a0058c125642e0064c817/a8cfefdd1ac031bbc125669b00493127/$FILE/MPI-I-98-1-022.ps |title=Fast Recursive Division |first=Christoph Burnikel |last=Joachim Ziegler |year=1998 |location=Max-Planck-Institut für Informatik |access-date=2021-09-10 |archive-date=2011-04-26 |archive-url=https://web.archive.org/web/20110426221250/http://domino.mpi-inf.mpg.de/internet/reports.nsf/efc044f1568a0058c125642e0064c817/a8cfefdd1ac031bbc125669b00493127/$FILE/MPI-I-98-1-022.ps |url-status=live }}</ref> や [[:en:Barrett reduction|Barrett reduction]] アルゴリズムがある<ref>{{Cite conference|url= http://portal.acm.org/citation.cfm?id=36688 |title=Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor |author=Paul Barrett |date=1987 |publisher=Springer-Verlag |book-title=Proceedings on Advances in cryptology---CRYPTO '86  |pages=311&ndash;323 |location=London, UK |isbn=0-387-18047-8}}</ref>。ニュートン法は同じ除数で複数の被除数に対して除算を行う場合に特に効率的で、除数の逆数を1度計算しておくと、毎回それを流用できる。

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

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

これを整数演算でおこなう場合は、256/3は当然正確な整数にはならないので、誤差があらわれる。しかし、シフト幅をより大きくし、値の範囲に注意すれば、常に不正確な部分はビットシフトによって捨てられる<ref>[http://gmplib.org/~tege/divcnst-pldi94.pdf Division by Invariant Integers using Multiplication] Torbjörn Granlund and Peter L. Montgomery. ACM SIGPLAN Notices Vol 29 Issue 6 (June 1994) 61–72</ref>ように変形できることがある。

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

場合によっては、定数による除算を一連のシフト操作と加減算に変換できることもある<ref>[http://techref.massmind.org/techref/method/math/divconst.htm Massmind: "Binary Division by a Constant"]</ref>。特に興味深いのは10による除算で、シフトと加減算で正確な商（と必要なら余り）が得られる<ref>R. A. Vowels, "Divide by 10", Australian Computer Journal (24), 1992, pp. 81-85.</ref>。

== 脚注 ==
{{Reflist}}

== 外部リンク ==
* [http://www.ecs.umass.edu/ece/koren/arith/simulator/ Computer Arithmetic Algorithms JavaScript Simulator] – 各種除算アルゴリズムのシミュレータがある。

{{DEFAULTSORT:しよさん}}
[[Category:コンピュータの算術]]
[[Category:除法]]