符号なし除算(回復法)
除算をFPGAに実装するとき、自分でRTLを書くのは難しそうなのでIP生成ツールを使う人が多いと思うが、最近はVitisHLSのような高位合成ツールがあるので、アルゴリズムを知っていると用途に合わせて除算器を作れる。そこで除算のアルゴリズムを合成可能なCコードで書いてみよう。ツールとしてVitisHLSの固定小数点型テンプレートap_int.hを使うが、SystemCのsystemc.hでも同様である。(→補足−固定小数点型テンプレートの使い方)
◇
まず、符号なし整数の除算を考える。なお固定小数点でも、最下位ビットを合わせれば整数の除算と同じである
|
a /d = q + e 分子a：符号なし、整数部IAビット 分母d：符号なし、整数部IDビット 商q： 符号なし、整数部IQビット、小数部FQビット (1) a,dが任意の値をとりうる場合はa/d≦2IA/1だからIQ=IAだけ必要 (2) a/d≦2IQ-2-FQとわかっている場合はビット幅IQ<IAで十分 例) N個の値xiの平均値はxiの値の範囲内に必ず収まる 誤差e：(1) 切捨の場合 0≦e＜2-FQ (2) 丸めの場合 -2-FQ-1≦e＜2-FQ-1 |
手順は、分子から分母を引いてみて、非負だったらq[i]=1にして次は引いた値を使い、負だったらq[i]=0にして次は引く前の値に戻せばよい。この方法は引き戻し法、回復法、Restoring法とよばれる。
簡単な例としてIA=3,ID=2,IQ=IA=3,FQ=0の場合で手計算すると次のようになる。
|
【回復法】 符号なし/符号なし |
説明 |
|
例）a[2:0]=101=(5) d[1:0]=11=(3) _______ 11) 101.0 r={00,a[2]} -) 11 d --- 110 t=r-d <0 q[2]=0 010 r={r[1:0],a[1]} -) 11 d --- 111 t=r-d <0 q[1]=0 101 r={r[1:0],a[0]} -) 11 d --- 010.0 t=r-d >=0 q[0]=1 10.0 r={t[1:0],0} -) 1.1 d -------- 00.1 t=r-d >=0 q[-1]=1 切捨の場合 q[2:0]=001 丸めの場合 q[2:0]=q[2:0]+q[-1]=010 |
←引き算はq[IA-1](aのMSB)から始める ←分母を引いて 負だったらq[i]=0にして次は引く前の値に戻す 非負だったらq[i]=1にして次は引いた値を使う t=r-d if(t<0) q[i]=0 r={r[1:0],a[i-1]} else q[i]=1 r={t[1:0],a[i-1]} 符号を判定するためt,rのビット幅は[ID:0]だけ必要
←切捨(0<=e<1)の場合はq[0]まで求めて終わり
←丸め(-0.5<=e<0.5)の場合はq[-1]まで求めて q[2:0]=q[2:0]+q[-1]とする |
◇◇
このアルゴリズムをVitisHLSで合成可能なCコードで書くと次のようになる。ここで商qを小数点以下まで求めるためにa*2FQ/d=qとaを左シフトしている。
|
#include <ap_int.h> /******************************************************** * division: a/d = q + e * param range: IQ<=IA * input range: 0<d,0<=a<=d*(2**(IQ)-2**(-FQ)) * error range: 0<=e<2**(-FQ) (truncation) *******************************************************/ template <int IA,int ID,int IQ,int FQ> ap_ufixed<IQ+FQ,IQ> div_uint_trn // (o) quotient (ap_uint<IA> a // (i) dividend ,ap_uint<ID> d // (i) divisor ) { ap_uint<IA+FQ> a_tmp=(ap_uint<IA+FQ>)a(IA-1,0)<<FQ; ap_uint<ID+1> r; ap_uint<ID+1> t; ap_ufixed<IQ+FQ,IQ> q; for(int i=IQ+FQ-1; i>=0; i--) { if(i==IQ+FQ-1) { r = a_tmp(IA+FQ-1,i); } else { r = (q[i+1]? t(ID-1,0):r(ID-1,0),a_tmp[i]); } t = r - d(ID-1,0); q[i] = !t[ID]; } return q; } |
|
#include <ap_int.h> /********************************************************* * division: a/d = q + e * param range: IQ<=IA * input range: 0<d,0<=a<=d*(2**(IQ)-2**(-FQ)) * error range: -2**(-FQ-1)<=e<2**(-FQ-1) (rounding) ********************************************************/ template <int IA,int ID,int IQ,int FQ> ap_ufixed<IQ+FQ,IQ> div_uint_rnd // (o) quotient (ap_uint<IA> a // (i) dividend ,ap_uint<ID> d // (i) divisor ) { ap_uint<IA+FQ+1> a_tmp=(ap_uint<IA+FQ+1>)a(IA-1,0)<<(FQ+1); ap_uint<ID+1> r; ap_uint<ID+1> t; ap_ufixed<IQ+FQ,IQ> q; for(int i=IQ+FQ; i>=0; i--) { if(i==IQ+FQ) { r = a_tmp(IA+FQ,i); } else { r = (q[i]? t(ID-1,0):r(ID-1,0),a_tmp[i]); } t = r - d(ID-1,0); if(i==0) { q(IQ+FQ-1,0) = q(IQ+FQ-1,0) + !t[ID]; } else { q[i-1] = !t[ID]; } } return q; } |
このコードは次のテストベンチで検証してある。
|
int main(int argc, char *argv[]) { const int IA=3; const int ID=3; #ifdef FULL_IQ const int IQ=3,FQ=1; // IQ=IA #else const int IQ=1,FQ=1; // IQ<IA #endif ap_uint<IA> a; ap_uint<ID> d; ap_ufixed<IQ+FQ,IQ> q; double a_val,d_val,q_val,e_val; printf("-- IA=%-d ID=%-d IQ=%-d FQ=%-d\n",IA,ID,IQ,FQ); /* a,dを全範囲振って商と誤差を表示する */ for(int j=1; j<pow(2,ID); j++) { for(int i=0; i<pow(2,IA); i++) { a = i; d = j; a_val = a.to_double(); d_val = d.to_double(); printf("-- a=%f d=%f",a_val,d_val); /* IQ<IA+1のとき大きいaはスキップする */ if(a_val<=d_val*(pow(2,IQ)-pow(2,-FQ))) { /* 切捨の場合 */ q = div_uint_trn<IA,ID,IQ,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]"); /* 丸めの場合 */ q = div_uint_rnd<IA,ID,IQ,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; } |