\begin{equation} y = f(x) = \frac{1}{\sqrt{x}} \end{equation} を高速に計算する。

float q_rsqrt(float number)
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;
 
  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;                       // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
 
  return y;
}

ニュートン法の1次近似であれば0x5f3759dfがマジックナンバーとなる。
2次近似だと変わったりするのかな

Cortex-M4での処理時間を計測してみる
ハードウェアのFPUを有効にする。浮動小数点は単精度浮動小数点(float 32bit)とする。
SystemCoreClockは170MHz、GCCの最適化オプションは-Og

// マジックナンバーを使っている有名な関数
float fast_inv_sqrt(float number){
 
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
 
    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
    // y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed
 
    return y;
}
 
// Cortex-M4コアで実装可能なARMの平方根計算用命令(vsqrt)
inline float cm4_mnemonic_sqrtf(float x){
  __asm__ ("vsqrt.f32 %0, %1" : "=t"(x) : "t"(x));
	return x;
}
 
void performanc_test_inv_sqrt(void){
  printf("i, time[ns], 1/sqrtf(), time[ns], fast_inv_sqrt, time[ns], 1/vsqrt()\r\n");
  uint32_t temp1, temp2, temp3;
  float temp1f, temp2f, temp3f;
  for(int16_t i = 1; i < 100; i++){
 
    reset_dwtcnt();
    temp1f = 1.0f / sqrtf((float)i);
    temp1 = get_dwtcnt(); // reset_dwtcnt()を実行してからの経過時間を記録
 
    reset_dwtcnt();
    temp2f = fast_inv_sqrt((float)i);
    temp2 = get_dwtcnt();
 
    reset_dwtcnt();
    temp3f = 1.0f / cm4_mnemonic_sqrtf((float)i);
    temp3 = get_dwtcnt();
 
    printf("%d, %ld,  %f, %ld, %f, %ld, %f\r\n", i, \
    calc_cycle_time_nanos_sec(temp1, SystemCoreClock), temp1f, \
    calc_cycle_time_nanos_sec(temp2, SystemCoreClock), temp2f, \
    calc_cycle_time_nanos_sec(temp3, SystemCoreClock), temp3f \
    );
  }
}

• 1/sqrtf() 558ns
• fast_inv_sqrt 335ns
• 1/vsqrt() 100ns

ARMコアに実装されている命令(VSQRT)を使ったほうが、1を割るための除算があったとしても有名な関数と比べると処理時間が短いことがわかった。

Fast inverse square root
Understanding Quake’s Fast Inverse Square Root
高速逆平方根(fast inverse square root)のアルゴリズム解説
FAST INVERSE SQUARE ROOT CHRIS LOMONT