開平(回復法)

[] [] []


平方根の手計算(開平)は学校で習わないが、知っておくとFPGAに実装できるので、開平のアルゴリズムを合成可能なCコードで書いてみよう

平方根は負数を考えなくてよいので、符号なし整数の開平を考える。一般の固定小数点でも、小数点の位置を移動すれば同じである。

sqrt(a) = q + e

元の数a:符号なし、整数部IAビット (0≦a≦2IA-1)

平方根q:符号なし、整数部IQビット、小数部FQビット (0≦q≦2IQ-2-FQ)

誤差e(1)切捨の場合 0≦e2-FQIQ=(IA+1)/2ビットだけ必要

     (a=(q+e)2(2IQ-2-FQ+2-FQ)2=22*IQ)

     (2)丸めの場合 -2-FQ-1≦e2-FQ-1 IQ=IA/2+1だけ必要

     (a=(q+e)2(2IQ-2-FQ+2-FQ-1)2=22*IQ-2IQ-FQ+2-2*FQ-2)

aの平方根qを求めることは

 a = q2 + r

となるqと余りrを求めることだから除算と似ている。ただ、除算では分母を引けばよいのに対して、開平では引く数をaから決めなければならない

今、整数qをビットごとにq=Σqj (qj=q[j]*2j ,q[j]=0,1)と表すと、恒等式

 (q3+q2+q1+q0)2 = q3q3 + (2q3+q2)q2 + (2(q3+q2)+q1)q1 + (2(q3+q2+q1)+q0)q0 …(1)

の左辺は大きな正方形、右辺の各項は小さな正方形とカギ形の面積に対応している。

 

そこで、右辺の第1項から順に、aから引いていけばよいことがわかる。

【開平】 a = q2 +r

引く数の計算    aからの引き算

   q3                a

+) q3            -) q3q3

  ---              -------

  2q3+q2             a-q3q3

+)    q2         -) q2(2q3+q2)

  ------            -----------

  2(q3+q2)+q1         a-q3q3-q2(2q3+q2)

+)        q1     -) q1(2(q3+q2)+q1)

  ----------       ------------------

  2(q3+q2+q1)+q0        a-q3q3-q2(2q3+q2)-q1(2(q3+q2)+q1)

+)           q0     -) q0(2(q3+q2+q1)+q0)

  -------------     --------------------------------

  2(q3+q2+q1+q0)      a-q3q3-q2(2q3+q2)-q1(2(q3+q2)+q1)-q0(2(q3+q2+q1)+q0)

結果

 q= q3+q2+q1+q0

 r= a-q3q3-q2(2q3+q2)-q1(2(q3+q2)+q1)- q0(2(q3+q2+q1)+q0)

実際の手順は、とりあえずq[i]=1として引いてみて、非負だったらq[i]=1で確定して次は引いた値を使い、負だったらq[i]=0に直して次は引く前の値に戻せばよい。(回復法、Restoring)

引く数は、例えばq[1]を求めるときは

 2(q3+q2)+q1 = 2*(q[3]*23+q[2]*22)+q[1]*21 =(4*q[3:2]+q[1])*21

と、既に求めたq[3:2]4倍して1を足したものになる(とりあえずq[1]=1とするので)

◇◇

簡単な例としてIA=6,IQ=3,FQ=0の場合で手計算してみると、次のようになる。

【回復法】 34 = (4+0+1)2 + 9

説明

例)a[5:0]=100010=(34)

   1        √0 10 00 10.00 r[4:2]={0,a[5:4]}

+) 1        -)0 01          {0,01}

  --          ----

  10          0 01          t[4:2] >=0

  101 q[2]=1    01 00       r[4:1]={t[3:2],a[3:2]}

+)  1         -)01 01       {0,q[2],01}

  ---           -----

  110           11 11       t[4:1] <0

  1001 q[1]=0    1 00 10    r[4:0]={r[3:1],a[1:0]}

+)   1         -)0 10 01    {0,q[2:1],01}

  ----           -------

  1010           0 10 01    t[4:0] >=0

  10101 q[0]=1     10 01.00 r[4:-1]={t[3:0],00}

+)    1          -)01 01 01 {0,q[2:0],01}

  -----            --------

  10100            00 11 11 t[4:-1] >=0

  10101 q[-1]=1

切捨の場合 q[2:0]=101

丸めの場合 q[2:0]=q[2:0]+q[-1]=110

 

引き算はq[IQ-1](aの上位2ビット)から始める

とりあえず01を引いて

 非負だったらq[i]=1にして次は引いた値を使い

 負だったらq[i]=0にして次は引く前の値に戻す

 if(t<0) q[i]=0 r[IQ+1:i]={r[IQ:i+1],a[2i+1:2i]}

 else    q[i]=1 r[IQ+1:i]={t[IQ:i+1],a[2i+1:2i]}

 符号を判定するためt,rのビット幅は[IQ+1:0]だけ必要

 次に引く数は4*q[IQ-1:i]+1になる

 

 

 

 

切捨の場合はq[0]まで求めて終わり

 

 

 

丸めの場合はq[-1]まで求めて

 q[2:0]=q[2:0]+q[-1]とする

◇◇◇

このアルゴリズムをVitisHLSで合成可能なCコードで書くと次のようになる。ここで平方根qを小数点以下まで求めるためにsqrt(a*22*FQ)=qaを左シフトしている。

#include <ap_int.h>

/*************************************************************

 * square root: sqrt(a) = q + e

 *   param range: IQ=(IA+1)/2

 *   input range: 0<=a

 *   error range: 0<=e<2**(-FQ) (truncation)

 ************************************************************/

template <int IA,int IQ,int FQ>

ap_ufixed<IQ+FQ,IQ> sqrt_uint_trn // (o) square root

  (ap_uint<IA> a        // (i) argument

  ) {

  ap_uint<2*(IQ+FQ)> a_tmp

           = (ap_uint<2*(IQ+FQ)>) a(IA-1,0) << 2*FQ;

  ap_uint<IQ+FQ+2> r;

  ap_uint<IQ+FQ+2> t;

  ap_ufixed<IQ+FQ,IQ> q;

  for(int i=IQ+FQ-1; i>=0; i--) {

    if(i==IQ+FQ-1) {

      r(IQ+FQ+1,i) = a_tmp(2*i+1,2*i);

      t(IQ+FQ+1,i) = r(IQ+FQ+1,i) - (bool)1;

    } else {

      r(IQ+FQ+1,i)

       =(q[i+1]? t(IQ+FQ,i+1):r(IQ+FQ,i+1),a_tmp(2*i+1,2*i));

      t(IQ+FQ+1,i)

       = r(IQ+FQ+1,i) - (q(IQ+FQ-1,i+1),(bool)0,(bool)1);

    }

    q[i] = !t[IQ+FQ+1];

  }

  return q;

}

#include <ap_int.h>

/*************************************************************

 * square root: sqrt(a) = q + e

 *   param range: IQ=IA/2+1

 *   input range: 0<=a

 *   error range: -2**(-FQ-1)<=e<2**(-FQ-1)(rounding)

 ************************************************************/

template <int IA,int IQ,int FQ>

ap_ufixed<IQ+FQ,IQ> sqrt_uint_rnd // (o) square root

  (ap_uint<IA> a        // (i) argument

  ) {

  ap_uint<2*(IQ+FQ+1)> a_tmp

           = (ap_uint<2*(IQ+FQ+1)>)a(IA-1,0) << 2*(FQ+1);

  ap_uint<IQ+FQ+3> r;

  ap_uint<IQ+FQ+3> t;

  ap_ufixed<IQ+FQ,IQ> q;

  for(int i=IQ+FQ; i>=0; i--) {

    if(i==IQ+FQ) {

      r(IQ+FQ+2,i) = a_tmp(2*i+1,2*i);

      t(IQ+FQ+2,i) = r(IQ+FQ+2,i) - (bool)1;

    } else {

      r(IQ+FQ+2,i)

       =(q[i]? t(IQ+FQ+1,i+1):r(IQ+FQ+1,i+1),a_tmp(2*i+1,2*i));

      t(IQ+FQ+2,i)

       = r(IQ+FQ+2,i) - (q(IQ+FQ-1,i),(bool)0,(bool)1);

    }

    if(i==0) {

      q(IQ+FQ-1,0) = q(IQ+FQ-1,0) + !t[IQ+FQ+2];

    } else {

      q[i-1] = !t[IQ+FQ+2];

    }

  }

  return q;

}

このコードは次のテストベンチで検証してある。

int main(int argc, char *argv[]) {

  const int IA=4;

#ifdef TRN

  const int IQ=(IA+1)/2,FQ=1;   // 切捨の場合

#else

  const int IQ=IA/2+1,FQ=1;     // 丸めの場合

#endif

  ap_uint<IA>         a;

  ap_ufixed<IQ+FQ,IQ> q;

  double a_val,q_val,e_val;

  printf("-- IA=%-d IQ=%-d FQ=%-d\n",IA,IQ,FQ);

  /* aを全範囲振って平方根と誤差を表示する */

  for(int i=0; i<pow(2,IA); i++) {

    a = i;

    a_val = a.to_double();

    printf("-- a=%f",a_val);

#ifdef TRN

    /* 切捨の場合 */

    q = sqrt_uint_trn<IA,IQ,FQ>(a);

    q_val = q.to_double();

    e_val = sqrt(a_val) - q_val;

    printf(" q_trn=%f e=%f",q_val,e_val);

    if(!(0<=e_val && e_val<pow(2,-FQ))) printf(" [FAIL]");

#else

    /* 丸めの場合 */

    q = sqrt_uint_rnd<IA,IQ,FQ>(a);

    q_val = q.to_double();

    e_val = sqrt(a_val) - q_val;

    printf(" q_rnd=%f e=%f",q_val,e_val);

    if(!(-pow(2,-FQ-1)<=e_val && e_val<pow(2,-FQ-1))) printf(" [FAIL]");

#endif

    printf("\n");

  }

  return 0;

}

◇◇◇

除算に非回復法があったように、開平でも非回復法が考えられる。恒等式

(q3+q2+q1+q0)2 = q3q3 + (2q3+q2)q2 + (2(q3+q2)+q1)q1 + (2(q3+q2+q1)+q0)q0 …(1)

において、各ビットをq[j]=0,1の代わりにq[j]=+1,-1とすればいいだけであるが、こうすると引く数の算出が面倒になるので、たぶん回路としては回復法よりも動作速度が遅くなると思う。


[] [] []