開平(回復法)
平方根の手計算(開平)は学校で習わないが、知っておくとFPGAに実装できるので、開平のアルゴリズムを合成可能なCコードで書いてみよう。
◇
平方根は負数を考えなくてよいので、符号なし整数の開平を考える。一般の固定小数点でも、小数点の位置を移動すれば同じである。
|
sqrt(a) = q + e 元の数a：符号なし、整数部IAビット (0≦a≦2IA-1) 平方根q：符号なし、整数部IQビット、小数部FQビット (0≦q≦2IQ-2-FQ) 誤差e：(1)切捨の場合 0≦e＜2-FQはIQ=(IA+1)/2ビットだけ必要 (∵a=(q+e)2＜(2IQ-2-FQ+2-FQ)2=22*IQ) (2)丸めの場合 -2-FQ-1≦e＜2-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)=qとaを左シフトしている。
|
#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とすればいいだけであるが、こうすると引く数の算出が面倒になるので、たぶん回路としては回復法よりも動作速度が遅くなると思う。