除算アルゴリズム

Method for division with remainder

除算アルゴリズムとは、2つの整数NとD（それぞれ分子と分母）が与えられ、それらの商と／または剰余（ユークリッド除算の結果）を計算するアルゴリズムです。手​​作業で適用されるものもあれば、デジタル回路設計やソフトウェアで利用されるものもあります。

除算アルゴリズムは、主に低速除算と高速除算の2つのカテゴリに分類されます。低速除算アルゴリズムは、反復ごとに最終的な商の1桁を生成します。低速除算の例としては、復元除算、非実行復元除算、非復元除算、SRT除算などがあります。高速除算法は、最終的な商の近似値から開始し、反復ごとに最終的な商の2倍の桁を生成します。^[1]ニュートン・ラプソン法とゴールドシュミット法はこのカテゴリに分類されます。

これらのアルゴリズムのバリエーションにより、高速な乗算アルゴリズムを使用できます。その結果、大きな整数の場合、どの乗算アルゴリズムを使用しても、除算に必要な計算時間は、定数倍まで乗算に必要な時間と変わりません。

議論はフォームを参照する。 ${\displaystyle N/D=(Q,R)}$

• N =分子（被除数）
• D =分母（除数）

は入力であり、

• Q =商
• R =余り

出力です。

繰り返し減算による除算

最も単純な除算アルゴリズムは、歴史的にはユークリッドの『原論』第 7 巻、命題 1で紹介された最大公約数アルゴリズムに組み込まれており、減算と比較のみを使用して 2 つの正の整数から余りを求めます。

function divide_unsigned ( N , D ) if D = 0 then error ( DivisionByZero ) end R : = N Q : = 0 while R ≥ D do R : = R − D Q : = Q + 1 end return ( Q , R ) end                                  

商と余りが存在し、一意であることの証明（ユークリッド除算で説明）により、加算、減算、比較を使用して、負の数と正の数の両方に適用できる完全な除算アルゴリズムが得られます。

function divide ( N , D ) if D = 0 then error ( DivisionByZero ) end if D < 0 then ( Q , R ) : = divide ( N , − D ) return ( − Q , R ) end if N < 0 then ( Q , R ) : = divide ( − N , D ) if R = 0 then return ( − Q , 0 ) else -- 例: N = -7、D = 3 -- divide(-N、D) = divide(7、3) = (2、1) -- R ≠ 0 なので、 return (-2 - 1、3 - 1) = (-3、2) -- チェック: (-3)*3 + 2 = -7 return ( − Q − 1 、D − R ) end end -- この時点で、 N ≥ 0 かつ D > 0 return divide_unsigned ( N 、D ) end                                                               

この手順は常にR ≥ 0を生成します。非常に単純ですが、Ω(Q)ステップを要するため、長除算のような低速な除算アルゴリズムよりも指数的に遅くなります。Qが小さいことが分かっている場合（出力に敏感なアルゴリズムであるため）、この手順は有用であり、実行可能な仕様として使用できます。

代替実装

別の実装では、剰余を増分し、除数に達したときにそれをリセットします。

の場合、アルゴリズムはが成り立つように計算します。 ただし、 ${\displaystyle x,y\in \mathbb {N} _{0}}$ ${\displaystyle q,r\,}$ ${\displaystyle x=qy+r\,}$ ${\displaystyle 0\leq r<y}$

${\displaystyle \operatorname {div} (x,y)=\left\lfloor x/y\right\rfloor }$

${\displaystyle x{\bmod {y}}=x-y\left\lfloor x/y\right\rfloor }$

Pythonの次のコードを考えてみます:

def divide_unsigned2 (分子: int ,分母: int ) -> tuple [ int , int ]商: int = 0剰余: int = 0 for _ in range (分子):剰余+= 1 if剰余==分母:商+= 1剰余= 0 return商、剰余                                   

注記

• 特別なケースは

${\displaystyle x{\bmod {1}}=0}$そして[ ^2]
${\displaystyle x{\bmod {0}}=x}$

• 変数quotientは決して読み込まれません。そのため、その代入（ハイライト表示）がコードから削除され、quotient出力リストからも削除された場合でも、divide_unsigned2 はdivide_unsigned と同様にを計算し続けます。 ${\displaystyle x{\bmod {y}}}$

カウンタマシンとしての代替実装
代替実装に基づく 単純なカウンタマシン（CM）が考えられます。CM命令は^{[3] [4]です。}

• Z (n): r _{n を}0 に置き換えます。
• S (n): r _nに1を加算します。
• J (m, n, q): r _m = r _nの場合はq 番目の命令にジャンプします。それ以外の場合は、プログラム内の次の命令に進みます。

カウンターマシンプログラムは^[5]

1: J(1,5,0)2: S(4)3: J(4,2,6)4: S(5)5: J(0,0,1)6: S(3)7: Z(4)8: S(5)9: J(0,0,1)

CMが初期レジスタ値R1=N、R2=D（残りのレジスタ=0）の計算を終えた後、

レジスタR3は、除算N/Dの商（整数部）を保持し、

レジスタ R4 には余りが保持されます。

長除算

長除法は、10進数で表された多桁の数を紙とペンで割る際に用いられる標準的なアルゴリズムです。このアルゴリズムでは、被除数の左端から右端へと段階的に移動させ、各段階で除数の（桁レベルでの）最大の倍数を減算します。そして、その倍数が商の桁となり、最終的な差が余りとなります。

この方法は、2進基数を用いると、以下の剰余付き（符号なし）整数除算アルゴリズムの基礎となります。短除算は、1桁の除数に適した長除算の短縮形です。チャンキング （部分商法またはハングマン法とも呼ばれる）は、長除算の効率は劣るものの、より理解しやすい方法です。各段階で現在持っている倍数よりも多くの倍数を減算できるようにすることで、より自由な形式の長除算の変種を開発することもできます。

剰余付き整数除算（符号なし）

以下のアルゴリズムは、有名な長除法の2進法バージョンであり、 N をDで割り、商をQに、余りをRに格納します。以下の擬似コードでは、すべての値は符号なし整数として扱われます。

D = 0の場合、error ( DivisionByZeroException )が発生します。Q : = 0 --商と余りをゼロに初期化します。R : = 0です。i : = n − 1 .. 0の場合、 doです。n は N のビット数です。R : = R << 1 -- R を 1 ビット左にシフトします。R ( 0 ) : = N ( i ) -- R の最下位ビットを分子のビット i と同じに設定します。R ≥ Dの場合、R : = R − D Q ( i ) : = 1です。end end                                             

例

_{N=1100 2} (12 ₁₀ )、D=100 ₂ (4 ₁₀ )とすると

ステップ1 : R=0、Q=0とする
ステップ2 : i=3とする（Nのビット数より1小さい）
ステップ3 : R=00とする（1だけ左シフトする）
ステップ4 : R=01とする（R(0)をN(i)とする）
ステップ5 : R < Dなので、ステートメントをスキップする

ステップ2 : i=2 に設定
ステップ3 : R=010
ステップ4 : R=011
ステップ5 : R < D、ステートメントはスキップ

ステップ2：i=1に設定
ステップ3：R=0110
ステップ4：R=0110
ステップ5：R>=D、ステートメントを入力
ステップ5b：R=10（R−D）
ステップ5c：Q=10（Q(i)を1に設定）

ステップ2：i=0に設定
ステップ3：R=100
ステップ4：R=100
ステップ5：R>=D、ステートメントを入力
ステップ5b：R=0（R−D）
ステップ5c：Q=11（Q(i)を1に設定）

終了
Q=11 ₂ (3 ₁₀ )およびR=0。

遅い除算法

低速除算法はすべて標準的な再帰方程式に基づいている^[6]

${\displaystyle R_{j+1}=B\times R_{j}-q_{n-(j+1)}\times D,}$

どこ：

• R _jは除算のj番目の部分剰余です（R(0) := N(i) ステップが含まれます）
• Bは基数（基数、コンピュータや電卓の内部では通常2）
• q _{n − ( j + 1)は商の}n −( j + 1)の位の数字であり、位の数字は最下位の 0 から最上位のn −1まで番号が付けられる。
• nは商の桁数です
• Dは除数である

分割の回復

復元除算は固定小数点小数に対して行われ、0 < D < Nという仮定に依存します。^{[引用が必要]}

商の数字qは数字セット{0,1}から形成されます。

2進数（基数2）の復元除算の基本アルゴリズムは次のとおりです。

R : = N D : = D << n -- R と D は N と Q の 2 倍のワード幅が必要ですfor i : = n − 1 .. 0 do -- たとえば 32 ビットの場合は 31..0 R : = 2 * R − D -- シフトされた値からの減算を試みます (2 を掛けると 2 進数表現ではシフトになります) if R >= 0 then q ( i ) : = 1 -- 結果ビット 1 else q ( i ) : = 0 -- 結果ビット 0 R : = R + D -- 新しい部分剰余は (復元された) シフトされた値ですend end                                             -- ここで、N = 分子、D = 分母、n = ビット数、R = 部分剰余、q(i) = 商のビット番号 i

実行されない復元除算は、2R の値が保存されることを除いて復元除算と似ています。そのため、R < 0 の場合にD を再度加算する必要はありません。

非復元部門

非復元除算では、商の桁として{0, 1}ではなく{−1, 1}の桁集合を用いる。このアルゴリズムはより複雑であるが、ハードウェア実装時には商ビットごとに判定と加減算が1回だけという利点がある。減算後に復元ステップがないため^[7] 、演算回数を最大半分まで削減し、高速化を図ることができる。^[8] 非負数の2進（基数2）非復元除算の基本アルゴリズムは以下の通りである。^[要検証]

R : = N D : = D << n -- R と D には、N と Q の 2 倍のワード幅が必要です。for i = n − 1 .. 0 do -- たとえば、32 ビットの場合は 31..0 if R >= 0 then q ( i ) : = + 1 R : = 2 * R − D else q ( i ) : = − 1 R : = 2 * R + D end if end -- 注: N = 分子、D = 分母、n = ビット数、R = 部分剰余、q(i) = 商のビット数。                                             

このアルゴリズムに従うと、商は-1と+1の数字からなる非標準形式になります。最終的な商を得るには、この形式を2進数に変換する必要があります。例：

次の商を数字セット{0,1}に変換します。
始める：	${\displaystyle Q=111{\bar {1}}1{\bar {1}}1{\bar {1}}}$
1. 正の項を形成する:	${\displaystyle P=11101010\,}$
2. 否定項をマスクする: [注 1]	${\displaystyle M=00010101\,}$
3. 引き算: : ${\displaystyle P-M}$	${\displaystyle Q=11010101\,}$
2の補数を使用せず、1の補数を使用する符号付き2進表記。

の -1 桁が、一般的であるようにゼロ (0) として格納されている場合、 となり、計算は簡単です。元の に対して 1 の補数 (ビットごとの補数) を実行します。 ${\displaystyle Q}$ ${\displaystyle P}$ ${\displaystyle Q}$ ${\displaystyle M}$ ${\displaystyle Q}$

Q : = Q − bit . bnot ( Q ) -- Q の −1 桁が通常どおりゼロとして表される場合は適切です。     

最後に、このアルゴリズムで計算される商は常に奇数であり、R の剰余は -D ≤ R < D の範囲にあります。例えば、5 / 2 = 3 R −1 です。剰余を正の値に変換するには、Q を非標準形式から標準形式に変換した後、1 回の復元ステップを実行します。

R < 0の場合Q : = Q − 1 R : = R + D -- 剰余が重要な場合のみ必要。                

実際の剰余は R >> n です。(除算の復元と同様に、R の下位ビットは商 Q のビットが生成されるのと同じ速度で消費されるため、両方に 1 つのシフト レジスタを使用するのが一般的です。)

SRT部門

SRT除算は、多くのマイクロプロセッサ実装において広く普及している除算手法である。^{[9] [10]このアルゴリズムは、} IBMのDW Sweeney、イリノイ大学のJames E. Robertson 、そしてロンドン大学インペリアル・カレッジのKD Tocherにちなんで名付けられた。彼らはほぼ同時期にこのアルゴリズムを独立して開発した（それぞれ1957年2月、1958年9月、1958年1月に発表）。^{[11] [12] [13]}

SRT 除算は非復元除算に似ていますが、被除数と除数に基づく参照テーブルを使用して商の各桁を決定します。

最も大きな違いは、商に冗長表現が使用されることです。例えば、基数4のSRT除算を実装する場合、商の各桁は5つの可能性（-2、-1、0、+1、+2）から選択されます。このため、商の桁の選択は完璧である必要はなく、後続の商の桁でわずかな誤差を修正できます。（例えば、商の桁のペア（0、+2）と（1、-2）は、0 × 4 + 2 = 1 × 4 − 2なので、同じです。）この許容範囲により、全角減算を必要とせず、被除数と除数の最上位ビットのみを使用して商の桁を選択できます。この簡素化により、2よりも高い基数を使用できるようになります。

非復元除算と同様に、最後の手順は、最後の商ビットを解決するための最後の全角減算と、商の標準の 2 進形式への変換です。

初代Intel Pentiumプロセッサの悪名高い浮動小数点除算バグは、ルックアップテーブルのコードが誤っていたために発生しました。Pentiumプロセッサは2048セルのテーブルを使用していましたが、そのうち1066セルに値が格納されるはずでしたが、5セルの値が誤って省略されていました。^{[14] [15] [16]}

高速除算法

ニュートン・ラプソン分割

ニュートン・ラプソン法では、ニュートン法を使用して の逆数を求め、その逆数に を掛けて最終的な商を求めます。 ${\displaystyle D}$ ${\displaystyle N}$ ${\displaystyle Q}$

ニュートン・ラプソン除算の手順は次のとおりです。

1. 除数 の逆数の推定値を計算します。 ${\displaystyle X_{0}}$ ${\displaystyle 1/D}$ ${\displaystyle D}$
2. 逆数のより正確な推定値を逐次的に計算する。ここでニュートン・ラプソン法が用いられる。 ${\displaystyle X_{1},X_{2},\ldots ,X_{S}}$
3. 被除数に除数の逆数を掛けて商を計算します。 ${\displaystyle Q=NX_{S}}$

ニュートン法を用いて の逆数を求めるには、で零点を持つ関数を見つける必要がある。そのような関数として明らかなのは だが、ニュートン・ラプソン反復法は の逆数が分かっていなければ計算できないため、役に立たない（しかも、反復的な改善ではなく、1ステップで正確な逆数を計算しようとする）。 は有効な関数であり、ニュートン・ラプソン反復法では以下の式が得られる。 ${\displaystyle D}$ ${\displaystyle f(X)}$ ${\displaystyle X=1/D}$ ${\displaystyle f(X)=DX-1}$ ${\displaystyle D}$ ${\displaystyle f(X)=(1/X)-D}$

$${\displaystyle 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}),}$$

これは、乗算と減算のみを使用して計算することも、2 つの融合した乗算と加算を使用して計算することもできます。 ${\displaystyle X_{i}}$

計算の観点から見ると、式と式は同等ではありません。2番目の式を用いて2nビットの精度の結果を得るには、と式の積を（nビット）の2倍の精度で計算する必要があります。 ^[要出典]一方、と式の積は、の先頭のnビット（2進小数点の後ろ）がゼロであるため、nビットの精度で計算するだけで済みます。 ${\displaystyle X_{i+1}=X_{i}+X_{i}(1-DX_{i})}$ ${\displaystyle X_{i+1}=X_{i}(2-DX_{i})}$ ${\displaystyle X_{i}}$ ${\displaystyle (2-DX_{i})}$ ${\displaystyle X_{i}}$ ${\displaystyle X_{i}}$ ${\displaystyle (1-DX_{i})}$ ${\displaystyle (1-DX_{i})}$

エラーが と定義されている場合、次のようになります。 ${\displaystyle \varepsilon _{i}=1-DX_{i}}$

${\displaystyle {\begin{aligned}\varepsilon _{i+1}&=1-DX_{i+1}\\&=1-D(X_{i}(2-DX_{i}))\\&=1-2DX_{i}+D^{2}X_{i}^{2}\\&=(1-DX_{i})^{2}\\&={\varepsilon _{i}}^{2}.\\\end{aligned}}}$

各反復ステップにおける誤差の二乗（いわゆるニュートン・ラプソン法の二次収束）は、結果の正しい桁数が反復ごとにほぼ倍増するという効果をもたらします。この特性は、対象となる数値の桁数が多い場合（例えば、大きな整数領域の場合）に非常に重要になります。しかし、これはまた、特に初期推定値が適切に選択されていない場合、この方法の初期収束が比較的遅くなる可能性があることも意味します。 ${\displaystyle X_{0}}$

初期見積もり

初期推定値を選択するという部分問題においては、除数Dにビットシフトを適用して0.5 ≤ D ≤ 1となるようにスケーリングする のが便利です。 同じビットシフトを分子Nに適用することで、商が変化しないことが保証されます。境界範囲内であれば、単純な多項式近似を用いて初期推定値を求めることができます。 ${\displaystyle X_{0}}$

区間における最悪の絶対誤差が最小となる線形近似は次のようになります。 ${\displaystyle [0.5,1]}$

${\displaystyle X_{0}={48 \over 17}-{32 \over 17}D.}$

線形近似の係数は次のように決定される。誤差の絶対値は である。誤差の絶対値の最大値の最小値は、にチェビシェフの等振動定理を適用することによって決定される。 の極小値はのときに発生し、その解は となる。この極小値における関数は、端点における関数と反対の符号、すなわち となる必要がある。2つの未知数における2つの方程式は、唯一の解と を持ち、最大誤差は である。この近似を用いると、初期値の誤差の絶対値は より小さくなる。 ${\displaystyle T_{0}+T_{1}D}$ ${\displaystyle |\varepsilon _{0}|=|1-D(T_{0}+T_{1}D)|}$ ${\displaystyle F(D)=1-D(T_{0}+T_{1}D)}$ ${\displaystyle F(D)}$ ${\displaystyle F'(D)=0}$ ${\displaystyle D=-T_{0}/(2T_{1})}$ ${\displaystyle F(1/2)=F(1)=-F(-T_{0}/(2T_{1}))}$ ${\displaystyle T_{0}=48/17}$ ${\displaystyle T_{1}=-32/17}$ ${\displaystyle F(1)=1/17}$

${\displaystyle \vert \varepsilon _{0}\vert \leq {1 \over 17}\approx 0.059.}$

区間内の最も適合度の高い二次曲線は ${\displaystyle 1/D}$

${\displaystyle X:={\frac {140}{33}}-{\frac {64}{11}}D+{\frac {256}{99}}D^{2}.}$

これは、誤差が第一種チェビシェフ多項式の再スケールされた三次多項式と等しくなるように選択され、誤差の絶対値が1/99以下になります。この改善はニュートン・ラプソン反復法と同等であり、計算コストは​​1反復未満です。 ${\displaystyle \log _{2}(\log 99/\log 17)\approx 0.7}$

Remezアルゴリズムを用いて係数を計算することで、2次以上の多項式近似を生成することが可能です。トレードオフとして、初期推定にはより多くの計算サイクルが必要になりますが、ニュートン・ラプソン法の反復回数は少なくなると期待できます。

この方法では収束は正確に2次式となるため、初期誤差から反復すると、 ${\displaystyle \varepsilon _{0}}$ ${\displaystyle S}$

${\displaystyle P=-2^{S}\log _{2}\varepsilon _{0}-1=2^{S}\log _{2}(1/\varepsilon _{0})-1}$

2進数。一般的な値は次のとおりです。

逆数精度の2進数

${\displaystyle \varepsilon _{0}}$	反復
	0	1	2	3	4
${\displaystyle 1/17}$	3.09	7.17	15.35	31.70	64.40
${\displaystyle 1/99}$	5.63	12.26	25.52	52.03	105.07

IEEE単精度では、二次初期推定値と2回の反復計算で十分な精度が得られますが、倍精度では3回の反復計算では限界があります。倍精度形式と倍精度拡張形式の両方では、線形初期推定値と4回の反復計算で十分です。

擬似コード

以下は、 NとDの商をPの 2 進数桁の精度で計算します。

D を M × 2 ^eとして表現します。ただし、1 ≤ M < 2 (標準浮動小数点表現)
D' := D / 2 ^e+1 // 0.5 から 1 の間でスケールし、ビットシフト/指数減算で実行できます
N' := N / 2 ^e+1
X := 48/17 − 32/17 × D' // D と同じ精度で定数を事前計算します
repeat times ${\displaystyle \left\lceil \log _{2}{\frac {P+1}{\log _{2}17}}\right\rceil \,}$ // 固定 P に基づいて事前計算できます
X := X + X × (1 - D' × X)
end
return N' × X

たとえば、倍精度浮動小数点除算の場合、この方法では 10 回の乗算、9 回の加算、および 2 回のシフトが使用されます。

3次反復

誤差を 3 乗するために 3 回の乗算を使用する反復があります。

${\displaystyle \varepsilon _{i}=1-DX_{i}}$

${\displaystyle Y_{i}=X_{i}\varepsilon _{i}}$

${\displaystyle X_{i+1}=X_{i}+Y_{i}+Y_{i}\varepsilon _{i}.}$

Y i _ε i_項は新しいものです。

上記を展開すると、次のように書ける。 ${\displaystyle X_{i+1}}$

${\displaystyle {\begin{aligned}X_{i+1}&=X_{i}+X_{i}\varepsilon _{i}+X_{i}\varepsilon _{i}^{2}\\&=X_{i}+X_{i}(1-DX_{i})+X_{i}(1-DX_{i})^{2}\\&=3X_{i}-3DX_{i}^{2}+D^{2}X_{i}^{3},\end{aligned}}}$

その結果、誤差項

${\displaystyle {\begin{aligned}\varepsilon _{i+1}&=1-DX_{i+1}\\&=1-3DX_{i}+3D^{2}X_{i}^{2}-D^{3}X_{i}^{3}\\&=(1-DX_{i})^{3}\\&=\varepsilon _{i}^{3}.\end{aligned}}}$

これは二次反復法の3/2の計算量ですが、収束性は同等であるため、わずかに効率が高くなります。言い換えれば、この手法を2回反復すると、誤差は9乗になりますが、二次反復法を3回反復すると誤差は8乗にしかならないのと同じ計算コストで済みます。 ${\displaystyle \log 3/\log 2\approx 1.585}$

反復後の正しいビットの数は ${\displaystyle S}$

${\displaystyle P=-3^{S}\log _{2}\varepsilon _{0}-1=3^{S}\log _{2}(1/\varepsilon _{0})-1}$

2進数。一般的な値は次のとおりです。

相互精度のビット

${\displaystyle \varepsilon _{0}}$	反復
	0	1	2	3
${\displaystyle 1/17}$	3.09	11.26	35.79	109.36
${\displaystyle 1/99}$	5.63	18.89	58.66	177.99

二次初期推定値と2回の三次反復計算を組み合わせることで、IEEE倍精度計算に十分な精度が得られます。二次反復計算と三次反復計算を併用することも可能です。

少なくとも1回の二次反復を使用すると、誤差が正になること、つまり逆数が過小評価されることが保証されます。^{[17] ：370}これにより、正確に丸められた商が必要な場合、後続の丸め手順が簡素化されます。

初期化または反復のいずれかで高次多項式を使用すると、必要な追加の乗算をより多くの反復の実行に費やした方がよいため、パフォーマンスが低下します。^{[引用が必要]}

ゴールドシュミット部門

ゴールドシュミット除算^[18]（ロバート・エリオット・ゴールドシュミットにちなんで）^{[19]は、共通の因数}F _{i を}被除数と除数の両方に繰り返し掛け合わせる反復処理を用いて、除数が1に収束するように選ばれる。これにより、被除数は求められる商Qに収束する。

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

ゴールドシュミット除算の手順は次のとおりです。

1. 乗算係数F _iの推定値を生成します。
2. 被除数と除数をF _iで乗算します。
3. 除数が 1 に十分近い場合は被除数を返し、そうでない場合はステップ 1 に戻ります。

N / Dが0 <  D < 1となるようにスケールされていると仮定すると 、各Fi_はDに基づきます。

${\displaystyle F_{i+1}=2-D_{i}.}$

被除数と除数を係数で乗算すると次のようになります。

${\displaystyle {\frac {N_{i+1}}{D_{i+1}}}={\frac {N_{i}}{D_{i}}}{\frac {F_{i+1}}{F_{i+1}}}.}$

十分な回数kの反復処理を行った後。 ${\displaystyle Q=N_{k}}$

ゴールドシュミット法はAMD Athlon CPUおよびそれ以降のモデルで使用されています。^{[20] [21]}これはアンダーソン・アール・ゴールドシュミット・パワーズ（AEGP）アルゴリズムとしても知られ、さまざまなIBMプロセッサに実装されています。^{[22] [23]}ニュートン・ラプソン法の実装と同じ速度で収束しますが、ゴールドシュミット法の利点の1つは、分子と分母の乗算を並列に実行できることです。^[23]

二項定理

ゴールドシュミット法は、二項定理によって簡略化できる因子を用いて用いることができます。 ⁠ ⁠ が ${\displaystyle N/D}$となるように2のべき乗でスケールされていると仮定します。 および とします。これは次式をもたらします。 ${\displaystyle D\in \left({\tfrac {1}{2}},1\right]}$ ${\displaystyle D=1-x}$ ${\displaystyle F_{i}=1+x^{2^{i}}}$

$${\displaystyle {\frac {N}{1-x}}={\frac {N\cdot (1+x)}{1-x^{2}}}={\frac {N\cdot (1+x)\cdot (1+x^{2})}{1-x^{4}}}=\cdots =Q'={\frac {N'=N\cdot (1+x)\cdot (1+x^{2})\cdot \cdot \cdot (1+x^{2^{(n-1)}})}{D'=1-x^{2^{n}}\approx 1}}}$$。

nステップ後、分母は次のように丸められる。 ${\displaystyle \left(x\in \left[0,{\tfrac {1}{2}}\right)\right)}$ ${\displaystyle 1-x^{2^{n}}}$相対誤差1

${\displaystyle \varepsilon _{n}={\frac {Q'-N'}{Q'}}=x^{2^{n}}}$

これはのときに最大となり、2進数の最小精度が提供されます。 ${\displaystyle 2^{-2^{n}}}$ ${\displaystyle x={\tfrac {1}{2}}}$ ${\displaystyle 2^{n}}$

大きな整数メソッド

ハードウェア実装用に設計された手法は、一般に、小数点以下数千または数百万の桁数を持つ整数には対応していません。このような状況は、暗号学におけるモジュラー縮約などで頻繁に発生します。このような大きな整数に対しては、より効率的な除算アルゴリズムによって、問題を少数の乗算を使用するように変換します。この乗算は、 Karatsuba アルゴリズム、Toom–Cook 乗算、Schönhage–Strassen アルゴリズムなどの漸近的に効率的な乗算アルゴリズムを使用して実行できます。その結果、除算の計算量は、乗算の計算量と同じオーダーになります (乗法定数まで)。例としては、前述のニュートン法による乗算への縮約、 ^[24]や、わずかに高速な Burnikel-Ziegler 除算、^[25] Barrett 縮約、Montgomery 縮約アルゴリズムなどがあります。^{[26] [検証が必要]}ニュートン法は、同じ除数で何度も除算しなければならない場合に特に効率的です。これは、最初のニュートン逆変換の後に、各除算に必要な乗算は1回だけ（切り捨て）であるためです。

定数による除算

定数Dによる除算は、その逆数による乗算と等価です。分母が定数であるため、その逆数 (1/ D ) も定数です。したがって、コンパイル時に(1/ D )の値を一度計算し、実行時には除算N/Dではなく乗算N ·(1/ D ) を実行することが可能です。浮動小数点演算では (1/ D )の使用にほとんど問題はありませんが^{[a] 、}整数演算では逆数は常にゼロになります (| D | > 1 と仮定)。

(1/ D ) を必ずしも使用する必要はありません。(1/ D )に簡約される任意の値 ( X / Y )を使用できます。例えば、3 で割る場合、1/3、2/6、3/9、または 194/582 を使用できます。したがって、Y が2 のべき乗であれば、除算手順は高速な右ビットシフトに簡約されます。N / D を ( N · X )/ Y として計算する効果は、除算を乗算とシフトに置き換えることです。N · ( X / Y ) はゼロに評価されるため、括弧が重要であることに注意してください。

しかし、D自体が2のべき乗でない限り、上記の条件を満たすXとYは存在しない。幸いなことに、( N・X )/ Yは、( X / Y )が1/ Dと完全に等しくなくても、整数演算においてN / Dと全く同じ結果を与える。ただし、近似によって生じる誤差はシフト演算で破棄されるビットに収まる程度に「十分近い」値である。^{[27] [28] [29]}バレット縮約は、Yの値を2のべき乗とすることで、Yによる除算を単純な右シフトにする。^[b]

具体的な固定小数点演算の例として、32ビットの符号なし整数の場合、3による除算は、⁠2863311531/2 ³³⁠ 、2863311531（ 16進数0xAAAAAAAB）を乗算し、その後33ビット右シフトします。2863311531の値は次のように計算されます。2 ³³/3⁠、その後切り上げられます。同様に、10で割る場合は、3435973837 (0xCCCCCCCD) を掛け、2 ³⁵ (または35 右ビットシフト) で割ることで表すことができます。^{[31] : p230-234} OEISでは、乗算の定数列をA346495、右シフトの定数列をA346496として提供しています。

除数D が2 の累乗ではない一般的なxビットの符号なし整数の除算では、次の恒等式により、 およびを事前計算した後、除算が 2 つのxビットの加算/減算、1 つのxビットxビットの乗算 (結果の上位半分のみを使用)、および複数のシフトに変換されます。 ${\displaystyle k=x+\lceil \log _{2}{D}\rceil }$ ${\displaystyle a=\left\lceil {\frac {2^{k}}{D}}\right\rceil -2^{x}}$

${\displaystyle \left\lfloor {\frac {N}{D}}\right\rfloor =\left\lfloor {\frac {\left\lfloor {\frac {N-b}{2}}\right\rfloor +b}{2^{k-x-1}}}\right\rfloor {\text{ where }}b=\left\lfloor {\frac {Na}{2^{x}}}\right\rfloor }$

場合によっては、「定数による除算」を一連のシフトと加算または減算に変換することで、さらに短時間で定数による除算を実行できます。^[32]特に興味深いのは10による除算です。この除算では正確な商が得られ、必要に応じて余りも得られます。^[33]

丸め誤差

除算を実行する際、正確な商 と余りはコンピュータの精度の限界内に収まるように近似されます。除算アルゴリズムは次のように定義されます。 ${\displaystyle q}$ ${\displaystyle r}$

${\displaystyle [a=bq+r]}$

どこ。 ${\displaystyle 0\leq r<|b|}$

浮動小数点演算では、商は、剰余はと表され、丸め誤差とが発生します。 ${\displaystyle q}$ ${\displaystyle {\tilde {q}}}$ ${\displaystyle r}$ ${\displaystyle {\tilde {r}}}$ ${\displaystyle \epsilon _{q}}$ ${\displaystyle \epsilon _{q}}$ ${\displaystyle \epsilon _{r}}$

${\displaystyle [{\tilde {q}}=q+\epsilon _{q}][{\tilde {r}}=r+\epsilon _{r}]}$

この丸めによって小さな誤差が生じ、それが後続の計算を通じて伝播し、蓄積される可能性があります。このような誤差は、反復処理やほぼ等しい値の減算において特に顕著になり、有意差が失われると言われています。これらの誤差を軽減するために、ガード桁の使用や高精度演算などの手法が用いられます。^{[34] [35]}

参照

• ギャレー部門
• 乗算アルゴリズム
• Pentium FDIVのバグ

注記

1. 最適化によって生じる問題がいかに「小さい」かにかかわらず、この逆最適化は不正確であるため、最新のコンパイラでは通常「高速演算」フラグの背後に隠されています。
2. 現代のコンパイラは一般的にこの整数乗算とシフトの最適化を実行します。ただし、実行時にのみ知られる定数の場合は、プログラム自身が最適化を実装する必要があります。^[30]

参考文献

1. Rodeheffer, Thomas L. (2008-08-26). ソフトウェア整数分割(PDF) (技術レポート). Microsoft Research, Silicon Valley.
2. 比較: プロシージャdivide_unsignedでは、厳密にdenominator正である必要があります。
3. Shepherdson, JC ; Sturgis, HE (1963). 「再帰関数の計算可能性」. Journal of the Association for Computing Machinery . 10 (2): 217– 255. doi : 10.1145/321160.321170 .
4. Cutland, Nigel (1980). 計算可能性：再帰関数理論入門(PDF) .ケンブリッジ大学出版局. p. 12. ISBN 9780521223843. 2025年9月16日閲覧。
5. このコードは「URMシミュレータ」上で動作します。ロサンゼルスのオクシデンタル大学。無制限レジスタマシン（URM）シミュレータ（エミュレータ）――「仮想URM」とも呼ばれます。これは、ケンブリッジ出版のナイジェル・J・カトランド著『計算可能性、再帰関数理論入門』に記載されているURM仕様をモデルにしています。
6. Morris, James E.; Iniewski, Krzysztof (2017-11-22). ナノエレクトロニクスデバイス応用ハンドブック. CRC Press. ISBN 978-1-351-83197-0。
7. Shaw, Robert F. (1950). 「バイナリコンピュータにおける算術演算」 . Review of Scientific Instruments . 21 (8): 690. Bibcode :1950RScI...21..687S. doi :10.1063/1.1745692. ISSN  0034-6748. 2022年2月28日時点のオリジナルよりアーカイブ。 2022年2月28日閲覧。
8. Flynn. 「Stanford EE486 (Advanced Computer Arithmetic Division) – Chapter 5 Handout (Division)」(PDF) .スタンフォード大学. 2022年4月18日時点のオリジナルよりアーカイブ(PDF) . 2019年6月24日閲覧。
9. Harris, David L.; Oberman, Stuart F.; Horowitz, Mark A. (1998年9月9日). SRT部門：アーキテクチャ、モデル、実装(PDF) (技術レポート). スタンフォード大学. 2016年12月24日時点のオリジナルよりアーカイブ(PDF) 。 2016年12月23日閲覧。
10. McCann, Mark; Pippenger, Nicholas (2005). 「SRT Division Algorithms as Dynamical Systems」. SIAM Journal on Computing . 34 (6): 1279– 1301. CiteSeerX 10.1.1.72.6993 . doi :10.1137/S009753970444106X. hdl :2429/12179. 2022年8月24日時点のオリジナルよりアーカイブ。2022年8月24日閲覧。 
11. Cocke, John; Sweeney, DW (1957年2月11日), High speed arithmetic in a parallel device (Company Memo), IBM, p. 20, 2022年8月24日時点のオリジナルよりアーカイブ。 2022年8月24日閲覧。{{citation}}: CS1 maint: location missing publisher (link)
12. Robertson, James (1958-09-01). 「新しいクラスのデジタル除算法」. IRE Transactions on Electronic Computers . EC-7 (3). IEEE: 218– 222. doi :10.1109/TEC.1958.5222579. hdl : 2027/uiuo.ark:/13960/t0gt7529c .
13. Tocher, KD (1958-01-01). 「自動バイナリコンピュータのための乗算と除算のテクニック」 . The Quarterly Journal of Mechanics and Applied Mathematics . 11 (3): 364– 384. doi :10.1093/qjmam/11.3.364. 2022年8月24日時点のオリジナルよりアーカイブ。 2022年8月24日閲覧。
14. 「浮動小数点欠陥の統計分析」Intel Corporation、1994年。2013年10月23日時点のオリジナルよりアーカイブ。 2013年10月22日閲覧。
15. Oberman, Stuart F.; Flynn, Michael J. (1995年7月). 除算アルゴリズムと実装の分析(PDF) (技術レポート). スタンフォード大学. CSL-TR-95-675. 2017年5月17日時点のオリジナルよりアーカイブ(PDF) 。 2016年12月23日閲覧。
16. Shirriff, Ken (2024年12月28日). 「Intelの4億7500万ドルのエラー：Pentium部門のバグの背後にあるシリコン」Righto . 2024年12月30日閲覧。
17. Ercegovac, Miloš D.; Lang, Tomás (2004). 「第7章 逆数、除算、逆平方根、反復近似法による平方根」. Digital Arithmetic . Morgan Kaufmann. pp.  367– 395. ISBN 1-55860-798-6。
18. Goldschmidt, Robert E. (1964). 収束による除算の応用(PDF) (論文). 理学修士論文. MIT OCLC 34136725. 2015年12月10日時点のオリジナルより アーカイブ(PDF) . 2015年9月15日閲覧。
19. 「著者」. IBM Journal of Research and Development . 11 : 125–127 . 1967. doi :10.1147/rd.111.0125. 2018年7月18日時点のオリジナルよりアーカイブ。
20. Oberman, Stuart F. (1999). 「AMD-K7マイクロプロセッサにおける浮動小数点除算および平方根アルゴリズムと実装」(PDF) . Proceedings 14th IEEE Symposium on Computer Arithmetic (Cat. No.99CB36336) . pp.  106– 115. doi :10.1109/ARITH.1999.762835. ISBN 0-7695-0116-8. S2CID  12793819. 2015年11月29日時点のオリジナルよりアーカイブ（PDF） . 2015年9月15日閲覧。
21. Soderquist, Peter; Leeser, Miriam (1997年7～8月). 「除算と平方根：適切な実装の選択」. IEEE Micro . 17 (4): 56– 66. doi :10.1109/40.612224.
22. SF Anderson、JG Earle、R.E. Goldschmidt、DM Powers. IBM 360/370 モデル91：浮動小数点演算ユニット、IBM Journal of Research and Development、1997年1月
23. Guy, Even; Peter, Siedel; Ferguson, Warren (2005年2月1日). 「ゴールドシュミットの除算アルゴリズムのパラメトリック誤差解析」. Journal of Computer and System Sciences . 70 (1): 118– 139. doi : 10.1016/j.jcss.2004.08.004 .
24. Hasselström, Karl (2003). 「大きな整数の高速除算：アルゴリズムの比較」（PDF）（コンピュータサイエンス修士論文）. 英国王立工科大学. 2017年7月8日時点のオリジナル（PDF）からアーカイブ。 2017年7月8日閲覧。
25. Joachim Ziegler、Christoph Burnikel (1998)、高速再帰部門、Max-Planck-Institut für Informatik、オリジナルから 2011 年 4 月 26 日にアーカイブ、2021 年 9 月 10 日に取得{{citation}}: CS1 maint: location missing publisher (link)
26. バレット、ポール (1987). 「標準デジタル信号プロセッサへのRivest ShamirとAdlemanの公開鍵暗号化アルゴリズムの実装」. Proceedings on Advances in Cryptology---CRYPTO '86 . ロンドン、英国: Springer-Verlag. pp.  311– 323. ISBN 0-387-18047-8。
27. Granlund, Torbjörn; Montgomery, Peter L. (1994年6月). 「乗算を用いた不変整数による除算」(PDF) . SIGPLAN Notices . 29 (6): 61– 72. CiteSeerX 10.1.1.1.2556 . doi :10.1145/773473.178249. 2019年6月6日時点のオリジナルよりアーカイブ(PDF) . 2015年12月8日閲覧. 
28. Möller, Niels; Granlund, Torbjörn (2011年2月). 「不変整数による除算の改良」(PDF) . IEEE Transactions on Computers . 60 (2): 165– 175. Bibcode :2011ITCmp..60..165M. doi :10.1109/TC.2010.143. S2CID 13347152. 2015年12月22日時点のオリジナルより アーカイブ( PDF) . 2015年12月8日閲覧.
29. thankful_fish. 「Labor of Division (Episode III): Faster Unsigned Division by Constants」2022年1月8日アーカイブ、Wayback Machine . 2011年。
30. foolish_fish. 「libdivide、最適化された整数除算」。2021年11月23日時点のオリジナルよりアーカイブ。 2021年7月6日閲覧。
31. ウォーレン・ジュニア、ヘンリー・S. (2013). Hacker's Delight (第2版). Addison Wesley - Pearson Education, Inc. ISBN  978-0-321-84268-8。
32. LaBudde, Robert A.; Golovchenko, Nikolai; Newton, James; Parker, David; Massmind: 「定数による二進除算」 2022年1月9日アーカイブ、Wayback Machine
33. Vowels, RA (1992). 「10で割る」.オーストラリアコンピュータジャーナル. 24 (3): 81–85 .
34. L. Popyack, Jeffrey (2000年6月). 「丸め誤差」.ドレクセル大学.
35. 「9. 機械数値、丸め誤差、および誤差伝播」チャールストン大学. 2021年2月8日.

さらに読む

• Savard, John JG (2018) [2006]. 「高度な算術技術」. quadibloc . 2018年7月3日時点のオリジナルよりアーカイブ。 2018年7月16日閲覧。

Retrieved from "https://en.wikipedia.org/w/index.php?title=Division_algorithm&oldid=1318866045"