SRT除算について

SRT除算は、回復法や非回復法と同様、商を上の桁から求めていくタイプの除算アルゴリズムです。 この時、回復法や非回復法と違い、商を冗長表現で求める点が特徴です。 商を冗長表現で求めることは、各桁の商の選択に多少の自由度を生み、これが回路的なメリットにつながります。

冗長表現とは、一つの数を表す方法が複数あるような記数法のことです。 N進法では通常、各桁の数字はN種類です。 しかしここでもっと多くの種類の数字を使えるようにすると、同じ数を表す方法が複数生まれます。 例えば、二進法で各桁の数字を0, 1, 2, 3から選ぶとすると、十三は1101とも1013とも213とも表せます(ほかにもいろいろあります)。

図1: 長除法(割り算の筆算)における部分商と部分剰余

商を冗長表現で求めることのメリットを考える前に、まずは普通の、冗長でない表現で商を求める場合を考えます。 冗長でない表現で商を求めるアルゴリズムである、回復法や非回復法を例にとってみます。 回復法や非回復法では、その時点での余り(部分剰余)を用いて各桁の商(部分商)を上の桁から順に決定していきます。 この時、部分商の決定を間違えることは許されません。 したがって、以下のように部分商を決定する際に部分剰余の正確な値が必要となります*1

  • 回復法では、除数と部分剰余を比べ、除数≧部分剰余なら部分商として1を立て、除数<部分剰余なら部分商として0を立てる。
  • 非回復法では、部分剰余の符号ビットが0なら部分商として1を立て、部分剰余の符号ビットが1なら部分商として-1を立てる。

商を冗長表現で求める場合、部分商を決定する際のある種の“間違い”を許容できます。 例えば、最終的な答えが十三の場合を考えてみます。 冗長でない二進法で求める場合、部分商の決定に間違いは許されず、順に1, 1, 0, 1と選んで1101とする必要があります。 一方、各桁に0, 1, 2, 3を使える冗長な二進法で求める場合、1, 0, と選んでしまっても許容範囲内です。 なぜなら、その後に2, 1と選んで1021としたり、あるいは1, 3と選んで1013としても、十三になるからです。 このように、特定の桁で多少“間違った”部分商を立ててしまっても、後続の部分商の選択をうまくやればリカバリーが可能です。

部分商を決定する際のある種の“間違い”が許容できることは、以下のような回路的メリットにつながります。

  • 部分商を決定する際、除数と部分剰余を全桁見る必要はなく、除数の上位数ビットと部分剰余の上位数ビットのみを見ればよい
    • テーブル引き*2での実装が可能になる。部分商を二~三桁同時に求めることも可能(Radix-4、Radix-8)*3
    • テーブル引きでなく回復法のように引き算した結果の符号を使って部分商を選択する実装の場合でも、その引き算を誤差あり(上位ビットだけの計算)で行うことができるようになる(参考:[2])
  • 部分剰余を正確に求める必要はなく、その上位数ビットのみが±1程度の誤差で求まればよい
    • 部分剰余は桁上げ保存加算器(と上位数ビットを確定させる小さな桁上げ先見加算器)で求めればよい。キャリーの伝搬を待つ必要がなくなる

さて、どの程度“間違った”部分商を立ててしまってもリカバリーが可能なのでしょうか。 これは、後続の桁の部分商をなるべく大きくした場合となるべく小さくした場合を考えることで求めることができます。

(なるべく大きくした場合)先の記数法では、後続の部分商を33333...とするのが最も大きくなります。 自分より下の桁が33333...となる時、その重みは自分の桁に対して0.33333...になります。 この値は三に等しいので、「真の商より(自分の桁で)三まで小さい商を立ててしまってもリカバリ可能」ということになります。

(なるべく小さくした場合)先の記数法では、後続の部分商を00000...とするのが最も小さくなります。 自分より下の桁が00000...となる時、その重みは自分の桁に対して0.00000...になります。 この値は零に等しいので、「真の商より(自分の桁で)零まで大きい商を立ててしまってもリカバリ可能(=ちょっとでも大きい商を立ててしまうとリカバリ不能)」ということになります。

ただし、これらの猶予は、部分剰余を正確に求めない場合には狭くなるので注意が必要です(PentiumのFDIVバグなど。参考:[2])。

ここまでがSRT除算の基本になります。 以下では、Radix-4・最小冗長表現・テーブル引き・部分剰余を正確に求める場合のSRT除算について詳しく見ていきます。 Radix-4とは四進で(=二進で二桁ずつ)商を求めること、最小冗長表現とはk進法に対して各桁の数字がk+1種類あることを意味します。

Radix-4・最小冗長表現・テーブル引き・部分剰余を正確に求める場合のSRT除算

Radix-4・最小冗長表現のSRT除算では、各桁の数字として-2, -1, 0, 1, 2を使います。 部分剰余を正確に求める場合、リカバリ可能な範囲を考えてみると、

  • 下の桁を22222...とすると最大で、その重みは自分の桁に対して0.22222...=三分の二
  • 下の桁を(-2)(-2)(-2)(-2)(-2)...とすると最小で、その重みは自分の桁に対してマイナス三分の二

となるので、真の商からプラスマイナス三分の二の範囲内に位置する商を立てればリカバリー可能です。 真の商からプラスマイナス二分の一の範囲内で商を立てなければいけなかった非回復法と比較してみてください。

以下では、除数・部分剰余・商の関係を見ていきます。 先ほどのリカバリ可能な範囲の議論から、除数Dと部分剰余R*4について、以下の関係が成り立ちます。

  • (4/3)D ≦ Rである時に限って、商として2を立ててよい。
  • (1/3)D ≦ R ≦ (5/3)Dである時に限って、商として1を立ててよい。
  • (-2/3)D ≦ R ≦ (2/3)Dである時に限って、商として0を立ててよい。
  • (-5/3)D ≦ R ≦ (-1/3)Dである時に限って、商として-1を立ててよい。
  • R ≦ (-4/3)Dである時に限って、商として-2を立ててよい。

ここで注目したいのは、これらの範囲には重複があることです。 例えば、(4/3)D ≦ R ≦ (5/3)Dの範囲の場合、商として1を立てても2を立ててもよいことになります。 こういった部分では、“間違い”が許容されるわけです。 一方、(2/3)D < R < (4/3)Dといった範囲では、“間違い”は許容されず、商として必ず1を立てなければなりません。

この重複の範囲をうまく利用すると、テーブル引きで商を求めることが可能です。 これを考えるために、まずはPDプロット*5を作ります。 横軸に除数D、縦軸に部分剰余Rを取って、先ほどの「(4/3)D ≦ R」などの境界線を引くと、図2のようになります*6。 ただし、除数Dは1≦D<2の範囲のみを図示してあります。 これは、この範囲内に正規化できることを考慮しています。 また、上下対称であるため、上半分のみを示しています。

図2: 最小冗長のRadix-4 SRTにおけるPDプロット

テーブル引きに使う表を作るには、これを格子状に区切り、各マスを商に対応した一つの色で塗ります。 それをやったのが図3のようになります(網掛け部分は、塗る色にいくらかの任意性があります)。 選択される商が許容される範囲内であることを確認してみてください。 境界線は斜めですが、猶予の部分をうまく活用することで、格子であっても許容された範囲内の商を選択することができています。 また、図3では、横軸を1/8刻み、縦軸を1/2刻みで区切ってあります。 除数Dは1≦D<2、部分剰余Rは-16/3<-(8/3)D≦R≦(8/3)D<16/3の範囲内であるため、これは除数Dの小数部の上位3bitと部分剰余Rの上位6bitのみでテーブル引きできることを意味します。

図3: 最小冗長のRadix-4 SRTで使用するテーブル構成法

これを使った浮動小数点数除算器のサンプルコードは以下になります。

int srt_table( int32_t rem, int32_t div ) {
  int d = ( div >> 20 ) - 8;
  assert( 0 <= d && d < 8 );

  int th12[] = { 6, 7, 8, 9, 9, 10, 11, 12 };
  int th01[] = { 2, 2, 3, 3, 3,  3,  4,  4 };

  rem >>= 21;
  if( rem < -th12[d] ) { return -2; }
  if( rem < -th01[d] ) { return -1; }
  if( rem <  th01[d] ) { return  0; }
  if( rem <  th12[d] ) { return  1; }
  return 2;
}

float srt_divider( float x, float y ) {
  assert( 1.0f <= x && x < 2.0f );
  assert( 1.0f <= y && y < 2.0f );

  int32_t X = x * 0x1.p+23;
  int32_t Y = y * 0x1.p+23;

  if( x < y ) { X *= 2; }

  int32_t quo = 0;
  int32_t rem = X;

  for( int i = 0; i < 13; ++i ) {
    const int q = srt_table( rem, Y );
    rem = ( rem - q * Y ) << 2;
    quo = ( quo << 2 ) + q;
  }
  // Here, quo has <1/3ULP error.
  
  // round nearest (never tie)
  if( quo & 1 ) {
    quo = rem < 0 ? quo - 1 : quo + 1;
  }

  float result = quo * ( x < y ? 0x1.p-25 : 0x1.p-24 );
  assert( result == x / y );
  return result;
}

参考文献

*1:非回復法では部分剰余の符号ビットだけ分かれば十分ですが、符号ビットを求めるのには全桁を求める場合とほとんど変わらない回路遅延が発生するので、「部分剰余の正確な値が必要」としました。

*2:テーブル引きと言っても、実用的にはROMを作るのではなく論理圧縮した回路とするようです。

*3:四桁(Intel PenrynのRadix-16 Dividerとか)や六桁(ARMのRadix-64 Floating-Point Dividerとか)を一サイクルで求めるものもありますが、これらは一サイクル内でRadix-4を複数回繰り返すという手法です。

*4:Partial reminderでPの記号を使う文献[1]もありますが、本稿ではより性質を表しているRを選択しました。文献[2]もRを使用しています。

*5:Partial reminderとdivisorでPDプロットというのが一般的[1][2]なようで、本稿でもそれに倣いました。文献[2]は部分剰余にRを使っているにもかかわらずPDプロットと呼んでいるので、RDプロットとは言わないのが普通のようです。

*6:この図は、PowerPointの6グリッドモードと8グリッドモードを駆使して作りました。