開平(回復法)

[] [] []


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

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

元の数a: 符号なし、整数部2Wビット

平方根q 符号なし、整数部Wビット、小数部Fビット (sqrt(22W)=2Wだからqの整数部a1/2ビット幅が必要)

誤差e = sqrt(a)-q  切捨て0 e 2-F の場合、aの範囲は任意

               丸め -2-F-1 e 2-F-1 の場合、aの範囲は a22W-2W-F

                (qmax=2W-2-Fだからa=(q+e)2(2W-2-F +2-F-1)2 = 22W-2W-F+2-2F-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[j]=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とするので)

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

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

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

                            10  手順

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

+) 1        -)0 01          ( 1)  {0,2'b01}

  --          ----

  10          0 01          ( 1)  t[2:0] >0

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

+)  1         -)01 01       ( 5)  {0,q[2],2'b01}

  ---           -----

  110           11 11       (-1)  t[3:0] <0

  1001 q[1]=0    1 00 10    (18)  r[4:0]={r[2:0],a[1:0]}

+)   1         -)0 10 01    ( 9)  {0,q[2:1],2'b01}

  ----           -------

  1010           0 10 01    ( 9)  t[4:0] >0

  10101 q[0]=1     10 01.00 (36)  r[5:0]={t[3:0],2'b00}

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

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

  10100            00 11 11 (15) t[5:0] >0

  10101 q[-1]=1

切捨の場合 q[2:0]=3'b101

丸めの場合 q[2:0]=q[2:0]+q[-1]=3'b110

 

 

 

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

 

 

 

if(t<0) q[i]=0 r[W+2-i:0]={r[W-i:0],a[i+1:i]} 引く前の値に戻す

 else    q[i]=1 r[W+2-i:0]={t[W-i:0],a[i+1:i]} 引いた値を使う

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

 

 

 

 

 

0<=r<2q+1(切捨)にするにはq[0]で終わり

 

 

 

-0.5<=r<0.5(丸め)にするにはq[-1]まで求めて

  さらにq[2:0]=q[2:0]+q[-1]とすればよい

◇◇

このアルゴリズムをVivadoHLS合成可能なCコードで書くと次のようになる。ここで平方根qを小数点以下まで求めるためにsqrt(a) =q*2-F すなわちsqrt(a*22F)=qが整数になるようにしている。またaは整数型ap_uint< >としたが、固定小数点型ap_ufixed< >に置き換えてもそのまま使えるように、小数点の位置に依存しないa.range( )で書いてある。

#include <ap_int.h>

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

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

 *   input range: 0<=a

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

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

template <int W,int F>

ap_ufixed<W+F,W> sqrt_uint_trn(ap_uint<2*W> a) {

    ap_uint<2*(W+F)> b = (ap_uint<2*(W+F)>) a.range(2*W-1,0) << 2*F;

    ap_uint<W+F+2> r;

    ap_int<W+F+2>  t;

    ap_ufixed<W+F,W> q;

    for(int i=W+F-1; i>=0; i--) {

        if(i==W+F-1) {

            r.range(W+F+1-i,0) = b.range(2*i+1,2*i);

            t.range(W+F+1-i,0) = r - 1;

        } else {

            r.range(W+F+1-i,0) = ((q[i+1]?t.range(W+F-1-i,0):r.range(W+F-1-i,0))

                                 ,b.range(2*i+1,2*i));

            t.range(W+F+1-i,0) = r - (q.range(W+F-1,i+1),(bool)0,(bool)1);

        }

        q[i] = !t[W+F+1-i];

    }

    return q;

}

#include <ap_int.h>

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

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

 *   input range: 0<=a<=2**(2W)-2**(W-F)

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

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

template <int W,int F>

ap_ufixed<W+F,W> sqrt_uint_rnd(ap_uint<2*W> a) {

    ap_uint<2*(W+F+1)> b = (ap_uint<2*(W+F+1)>) a.range(2*W-1,0) << 2*(F+1);

    ap_uint<W+F+3> r;

    ap_int<W+F+3>  t;

    ap_ufixed<W+F,W> q;

    for(int i=W+F; i>=0; i--) {

        if(i==W+F) {

            r.range(W+F+2-i,0) = b(2*i+1,2*i);

            t.range(W+F+2-i,0) = r - 1;

        } else {

            r.range(W+F+2-i,0) = ((q[i]? t.range(W+F-i,0):r.range(W+F-i,0))

                                 ,b(2*i+1,2*i));

            t.range(W+F+2-i,0) = r - (q.range(W+F-1,i),(bool)0,(bool)1);

        }

        if(i==0) {

            q.range(W+F-1,0) = q.range(W+F-1,0) + !t[W+F+2];

        } else {

            q[i-1] = !t[W+F+2-i];

        }

    }

    return q;

}

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

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

    const int W = 2;

    const int F = 0;

    ap_uint<2*W>     a;

    ap_ufixed<W+F,W> q;

    double e;

   

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

        a = i;

        cout << "-- a=" << setw(3) << a;

    /* 切捨ての場合 */

        q = sqrt_uint_trn<W,F>(a);

        e = sqrt(a.to_double()) - q.to_double();

        if(!(0<=e && e<pow(2,-F))) cout << " FAIL";

        cout << " q_trn=" << setw(3) << q << " e=" << setw(10) << e;

        /* 丸めの場合 */

        q = sqrt_uint_rnd<W,F>(a);

        e = sqrt(a.to_double()) - q.to_double();

        cout << " q_rnd=" << setw(3) << q << " e=" << setw(10) << e;

        if(!(-pow(2,-F-1)<=e && e<pow(2,-F-1))) cout << " FAIL";

        cout << endl;

    }

    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とすればいいだけであるが、こうすると引く数の算出が面倒になるので、たぶん回路としては回復法よりも動作速度が遅くなると思う。


[] [] []