正規化の除算(非回復法)

[] [] []


ベクトルの内積を正規化するときは、商が±1以内に限られるので、通常より少ない引き算(小さい回路)で除算ができる。

-1 ≦ ab / |a||b|=Σaj bj / sqrt(Σa2Σb2 ) +1

すなわち通常は被除数のMSBから始める引き算を1の位から始めればよい。また除数が正だから、非回復法で被除数の正負のみによって引き算か足し算かを決められる。

実際に正規化するときは、除数dは平方根sqrt(Σa2Σb2 )で計算するため誤差がでて-1q+1の範囲を少しオーバーする可能性があるが、そうなっても大丈夫である。なぜなら、商はq[1]まで求まるので-2q2-2-Fの範囲で使えるからである。

被除数a: 符号付き、整数部W+1ビット

除数d:  符号なし、整数部Wビット

q    符号付き、整数部2ビット、小数部Fビット。(-1q+1の範囲まで求めるためqの整数部2ビット必要)

剰余r = a-d*q  切捨て0 rd*2-F の場合、a,dの範囲は-2da+2d

             丸め -d*2-F-1 rd*2-F-1の場合、a,dの範囲は-2da(2-2-F-1)*d

              (∵ -2q2-2-Fだからa=d*q+r(2-2-F+2-F-1)*d=(2-2-F-1)*d )

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

【非回復法】負/

例)a[2:0]=101 (-3)  d[1:0]=11 ( 3)

          _______  10  手順

       11)1101.0

          1101     (-3)  r={符号拡張a[2:0]}

q[0]=-1 +)0011     ( 3)  d

          ----

          0000     ( 0)  t=r+d

           000.0   ( 0)  r={t[1:0],0}

q[-1]=+1 -)001.1   ( 3)  d

           -----

           110.1   (-3)  t=r-d

            10.10  (-6)  r={t[1:0],0}

q[-2]=-1  +)001.1  ( 3)  d

            -----

            110.1  (-3)  t=r-d

 

切捨の場合q[0:-1]=-1+0.5=-0.5-1.0

丸めの場合q[0:-2]=-1+0.5-0.25=-0.75-1.0

 

 

 

←引き算はq[0]から始める

if(r<0) q[i]=-1 t=r+d  逆符号にして引く

 else    q[i]=+1 t=r-d  そのまま引く

 

 

 

0<=r<d*2-1(切捨て)にするにはq[-1]まで求めて終わり

 +1のビットをp[i]とおくと

 q[1:-1] = 3'b100+{p[0:-1],1'b1}

 t<0(引きすぎ)のときは、q=q-2-1とすればよいので

 if(t<0) q[1:-1]=3'b100+{p[0:-1],1'b0}

 else    q[1:-1]=3'b100+{p[0:-1],1'b1}

-d*2-2<=r<d*2-2 (丸め)にするにはq[-2]まで求めて

 if(t<0) q[1:-2]=4'b1000+{p[0:-2],1'b0}

 else    q[1:-2]=4'b1000+{p[0:-2],1'b1}

 q[1:-1]=q[1:-1]+q[-2](切上げ)とすればよいので

 if(t<0) q[1:-1]=3'b100+p[0:-2]

 else    q[1:-1]=3'b100+p[0:-2]+1'b1

◇◇

このアルゴリズムをVivadoHLS合成可能なCコードで書くと次のようになる。ここで、a,dは整数型ap_int< >としたが、固定小数点型ap_fixed< >に置き換えてもそのまま使えるように、小数点の位置に依存しないa.range( )d.range( )で書いてある。

#include <ap_int.h>

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

 * division a/d = q + e

 *   input range: -2d<=a<2d

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

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

template <int W,int F>

ap_fixed<F+2,2> div_norm_trn

    (ap_int<W+1> a      // (i) dividend

    ,ap_uint<W>  d      // (i) divisor

    ) {

    ap_int<W+2> r;

    ap_int<W+2> t;

    ap_uint<F+1> p;

    ap_fixed<F+2,2> q;  // (o) quotient

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

        if(i==F) {

            r = (ap_int<W+1>) a.range(W,0);

        } else {

            r = (t.range(W,0),(bool)0);

        }

        p[i] = !r[W+1];

        t = (r[W+1])? r + d.range(W-1,0) : r - d.range(W-1,0);

    }

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

    return q;

}

#include <ap_int.h>

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

 * division a/d = q + e

 *   input range: -2d<=a<(2-2**(-F-1))*d

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

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

template <int W,int F>

ap_fixed<F+2,2> div_norm_rnd

    (ap_int<W+1> a      // (i) dividend

    ,ap_uint<W>  d      // (i) divisor

    ) {

    ap_int<W+2> r;

    ap_int<W+2> t;

    ap_uint<F+2> p;

    ap_fixed<F+2,2> q;  // (o) quotient

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

        if(i==F+1) {

            r = (ap_int<W+1>) a.range(W,0);

        } else {

            r = (t.range(W,0),(bool)0);

        }

        p[i] = !r[W+1];

        t = (r[W+1])? r + d.range(W-1,0) : r - d.range(W-1,0);

    }

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

    return q;

}

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

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

    const int W = 3;

    const int F = 1;

    ap_int<W+1> a;

    ap_uint<W>  d;

    ap_fixed<F+2,2> q;

    double e;

    int i_min,i_max;

    for(int j=1; j<pow(2,W); j++) {

        i_min = (-pow(2,W)<-2*j)? -2*j : -pow(2,W);

        i_max = ( pow(2,W)> 2*j)?  2*j :  pow(2,W);

        for(int i=i_min; i<i_max; i++) {

            a = i;

            d = j;

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

      /* 切捨ての場合 */

            q = div_norm_trn<W,F>(a,d);

            e = a.to_double() / d.to_double() - q.to_double();

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

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

            /* 丸めの場合 */

            q = div_norm_rnd<W,F>(a,d);

            e = a.to_double() / d.to_double() - q.to_double();

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

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

            cout << endl;

        }

    }

    return 0;

}


[] [] []