開平(回復法)

[←]　[↑]　[→]

---

符号なし整数の場合の平方根の手計算(開平)は開平(回復法)に書いたとおりで、これを固定小数点に直したものは次のとおり。誤差は切捨でなく丸め(四捨五入)とする。

sqrt(a) = q + e

元の数a：符号なし、整数部IAビット、小数部FAビット (0≦a≦2^IA-2^-FA)

　　　 　FQ=FA/2なのでFAは2の倍数が必要

平方根q：符号なし、整数部IQビット、小数部FQビット (0≦q≦2^IQ-2^-FQ)

誤差e：丸めの場合 -2^-FQ-1≦e＜2^-FQ-1は IQ=IA/2+1だけ必要

　　　　(∵a=(q+e)²＜(2^IQ-2^-FQ+2^-FQ-1)²=2^2*IQ-2^IQ-FQ+2^-2*FQ-2)

VitisHLSで合成可能なCコードは

#include <ap_int.h>

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

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

 *   param range: IQ=IA/2+1 FQ>=FA/2 FA=0,2,4,...

 *   input range: 0<=a

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

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

template <int IA,int FA,int IQ,int FQ>

ap_ufixed<IQ+FQ,IQ> sqrt_ufixed // (o) square root

  (ap_ufixed<IA+FA,IA> a        // (i) argument

  ) {

  ap_uint<2*(IQ+FQ+1)> a_tmp

    = (ap_uint<2*(IQ+FQ+1)>) a(IA+FA-1,0) << 2*(FQ-FA/2+1);

  ap_uint<IQ+FQ+3> r;

  ap_int<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,0) = a_tmp(2*i+1,2*i);

      t(IQ+FQ+2-i,0) = r - 1;

    } else {

      r(IQ+FQ+2-i,0)

        = (q[i]? t(IQ+FQ-i,0):r(IQ+FQ-i,0),a_tmp(2*i+1,2*i));

      t(IQ+FQ+2-i,0) = r - (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-i];

    }

  }

  return q;

}

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

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

  const int IA=4,FA=2;  // FA=0,2,4,...

#ifdef JUST_FQ

  const int IQ=IA/2+1,FQ=1;     // FQ=FA/2

#else

  const int IQ=IA/2+1,FQ=3;     // FQ>FA/2

#endif

  ap_ufixed<IA+FA,IA> a;

  ap_ufixed<IQ+FQ,IQ> q;

  double a_val,q_val,e_val;

  printf("-- IA=%-d FA=%-d IQ=%-d FQ=%-d\n",IA,FA,IQ,FQ);

/* aを全範囲振って平方根と誤差を表示する */

  for(int i=0; i<pow(2,IA+FA); i++) {

    a = (float) i * pow(2,-FA);

    a_val = a.to_double();

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

    q = sqrt_ufixed<IA,FA,IQ,FQ>(a);

    q_val = q.to_double();

    e_val = sqrt(a_val) - q_val;

    printf(" q=%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;

}

---

[←]　[↑]　[→]