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

[←]　[↑]　[→]

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

-1≦ a・b / |a||b|=Σa_jb_j / sqrt(Σa_j²Σb_j² ) ≦+1

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

実際に正規化するときは、分母dは平方根sqrt(Σa_j²Σb_j²)で計算するため誤差がでて-1≦q≦+1を少しはみだす可能性があるが、商は整数部2ビットまで求めている(-2≦q≦2-2^-FQ)ので大丈夫である。

a = d * q + e

分子a：符号つき、整数部ID+1ビット

分母d：符号なし、整数部IDビット

商q：　符号つき、整数部2ビット、小数部FQビット (-2≦q≦2-2^-FQ)

　　　　(∵-1≦q≦1の範囲まで求めるためqの整数部は2ビット必要)

誤差e：(1)切捨の場合 0 ≦e＜2^-FQ、a,dの範囲は-2≦a/d＜2

　　　　　(∵a/d=q+e＜2-2^-FQ+2^-FQ=2 )

　　　 (2)丸めの場合 -2^-FQ-1≦e＜2^-FQ-1、a,dの範囲は-2≦a/d＜(2-2^-FQ-1)

　　　　　(∵a/d=q+e＜2-2^-FQ+2^-FQ-1=2-2^-FQ-1 )

◇

簡単な例としてID=2,FQ=1の場合で手計算すると次のようになる。

【非回復法】負/正

説明

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

_______

11)1101.00 r={符号拡張a[2:0]}

q[0]=-1 +)0011 d

----

0000 t=r+d >=0

000.0 r={t[2:0],0}

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

-----

110.1 t=r-d <0

10.10 r={t[2:0],0}

q[-2]=-1 +)00.11 d

-----

11.01 t=r-d <0

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

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

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

←負だったらq[i]=-1として分母を足す

　非負だったらq[i]=+1として分母を引く

　if(r<0) q[i]=-1 t=r+d

　else　　q[i]=+1 t=r-d

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

←切捨(0<=e<0.5)の場合はq[-1]まで求めて終わり

　+1のビットをp[i]とおくとq[1:-1]=100+{p[0:-1],1}

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

　if(t<0) q[1:-1]=100+{p[0:-1],0}

　else q[1:-1]=100+{p[0:-1],1}

←丸め(-0.25<=e<0.25)の場合はq[-2]まで求めて

　if(t<0) q[1:-2]=1000+{p[0:-2],0}

　else q[1:-2]=1000+{p[0:-2],1}

　q[1:-1]=q[1:-1]+q[-2]とすればよいので

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

　else q[1:-1]=100+p[0:-2]+1

◇◇

このアルゴリズムをVitisHLSで合成可能なCコードで書くと次のようになる。

#include <ap_int.h>

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

* division: a/d = q + e

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

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

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

template <int ID,int FQ>

ap_fixed<2+FQ,2> div_norm_trn // (o) quotient

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

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

) {

ap_int<ID+2> r;

ap_int<ID+2> t;

ap_uint<1+FQ> p;

ap_fixed<2+FQ,2> q;

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

if(i==FQ) {

r = (ap_int<ID+1>) a(ID,0); // sign-extension

} else {

r = (t(ID,0),(bool)0);

}

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

t = (p[i])? r - d(ID-1,0) : r + d(ID-1,0);

}

q(1+FQ,0)=((ap_uint<2+FQ>)1<<(1+FQ))+(p,!t[ID+1]);

return q;

}

#include <ap_int.h>

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

* division: a/d = q + e

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

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

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

template <int ID,int FQ>

ap_fixed<2+FQ,2> div_norm_rnd // (o) quotient

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

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

) {

ap_int<ID+2> r;

ap_int<ID+2> t;

ap_uint<2+FQ> p;

ap_fixed<2+FQ,2> q;

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

if(i==1+FQ) {

r = (ap_int<ID+1>) a(ID,0); // sign-extension

} else {

r = (t(ID,0),(bool)0);

}

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

t = (p[i])? r - d(ID-1,0) : r + d(ID-1,0);

}

q(1+FQ,0)= ((ap_uint<2+FQ>)1<<(1+FQ))+ p+ !t[ID+1];

return q;

}

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

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

const int ID = 3;

const int FQ = 1;

ap_int<ID+1> a;

ap_uint<ID> d;

ap_fixed<2+FQ,2> q;

double a_val,d_val,q_val,e_val;

/* a,dを全範囲振って商と誤差を表示する */

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

for(int i=-2*j; i<2*j; i++) {

a = i;

d = j;

a_val = a.to_double();

d_val = d.to_double();

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

/* 切捨の場合 */

q = div_norm_trn<ID,FQ>(a,d);

q_val = q.to_double();

e_val = a_val/d_val - q_val;

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

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

/* 丸めの場合,入力範囲外のaはスキップする */

if(a_val < d_val*(2-pow(2,-FQ-1))) {

q = div_norm_rnd<ID,FQ>(a,d);

q_val = q.to_double();

e_val = a_val/d_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]");

}

printf("\n");

}

return 0;

}

[←]　[↑]　[→]