数値演算のアルゴリズム

2002.12.09 Yutaka Yoshisaka.

３Ｄの計算処理では、「正確性よりも速度を求める」という場合がよくあります（特にリアルタイム）。
そのあたりで使えそうな、数値演算のアルゴリズムをまとめてみました。

sqrt：平方根を求める

C言語では"math.h"の「sqrt」で平方根を計算します。
これと同等の機能をするプログラムは以下のようになります。

浮動小数点での平方根の計算

double mySqrt(double x)
{
  double s,last;

  if(x<=0.0) return 0.0;

  if (x > 1) s = x;  else s = 1;
  do {
    last = s;
    s = (x / s + s) * 0.5;
  } while (s < last);

  return last;
}

整数での平方根の計算

int myISqrt(int x)
{
  int s, t;

  if (x <= 0) return 0;

  s = 1;  t = x;
  while (s < t) {  s <<= 1;  t >>= 1;  }
  do {
    t = s;
    s = (x / s + s) >> 1;
  } while (s < t);

  return t;
}

Windows環境・Mac OS X環境では、すでに標準の「sqrt( )」のほうが速度が出てますので、出番はないでしょう(^_^;;。
DOS時代は、整数での平方根演算は上の方法のほうが速かったです。浮動小数点演算コプロセッサがついていなかった過去の話です。
しかし、浮動小数点演算が弱いPS2などでは、まだこの整数平方根演算のほうが速いみたいです(^_^;;

rsqrt：逆数平方根を求める

「1.0/sqrt(x)」の計算です。
これは、以下の方法で近似したほうが圧倒的に速いです。

float myRSqrt(float x)
{
  unsigned int iDat;
  double dx,dy;

  iDat = ((0x7f000000 + 0x3f800000) - (*((unsigned int *)&x))) >> 1;
  dy = *((float *)&iDat);

  dx = x;
  dy = dy * 1.47 - dx * 0.47 * dy * dy * dy;
  dy = (3.0 - dy * dy * dx) * 0.5 * dy;

  return (float)dy;
}

VC++では、コンパイラ時にケタオーバーのエラーが出るのですが無視しても結構です。
除算も平方根も使うことなく、「1.0/sqrt(x)」が計算できてしまいます。不思議ですね、私もよく分かりません(^_^;;
ただ、上記の方法では値によっては、計算時のオーバーフロー・アンダーフローの関係上、動作が挙動不審になったりする現象が起きていました。

そこで、以下のように補正するようにします。

float myRSqrt(float x)
{
  unsigned int iDat;
  double dx,dy;

  if(x<(1e-5)) return 0.0f;

  if(x>(1e+5)){

    //少数部分をカット
#if _WINDOWS
    x = (float)floor(x);
#else
    x = (float)((int)x);
#endif

    x = (float)((double)x * 0.00390625);  // x=x/256.0

    iDat = ((0x7f000000 + 0x3f800000) - (*((unsigned int *)&x))) >> 1;
    dy = *((float *)&iDat);

    dx = x;
    dy = dy * 1.47 - dx * 0.47 * dy * dy * dy;
    dy = (3.0 - dy * dy * dx) * 0.03125 * dy;

  } else {

    iDat = ((0x7f000000 + 0x3f800000) - (*((unsigned int *)&x))) >> 1;
    dy = *((float *)&iDat);

    dx = x;
    dy = dy * 1.47 - dx * 0.47 * dy * dy * dy;
    dy = (3.0 - dy * dy * dx) * 0.5 * dy;
  }

  return (float)dy;
}

「x=0.00001」以下の場合は、0とみなしてスキップしてます。
また、「x=100000」以上の場合はケタあふれが出ると想定して、先に256で割って値を小さくしてしまいます。
そして、逆平方根を計算後にsqrt(256)=16で割っています（0.5/16 ==> 0.03125）。
なお、アンダーフロー・オーバーフローが出るのを防ぐために、double計算している点に注意してください。

「#if _WINDOWS」で、Windows環境では「floor」にて整数部のみ取り出してます。
これは、最適化によっては「x = (float)((int)x);」が省略されてしまう可能性があるのと、 floor関数自身が、少ないクロック数で実行できるためにこのようにしています。
OS Xでは、「floor」自身が遅いため（また、最適化でも省略されないため「x = (float)((int)x);」としています。

ただ、注意点として「1.0/sqrt(x)よりは、誤差が大きい」ところに注意してください。
無視できる程度といえば、そのような気もします。

平方根の計算は、この逆数平方根を用いて「x * rsqrt(x)」として求めてしまうのも手ですね。
OS Xでは、普通の平方根よりも、この方が速いです。

hypot：(x,y)の長さを近似する

double myHypot(double x,double y)
{
  double len,len1,len2;

  x = (x < 0) ? -x : x;
  y = (y < 0) ? -y : y;
  if(x >= y){
    len1 = x;
    len2 = y;
  } else {
    len1 = y;
    len2 = x;
  }

  len = 0.9604 * len1 + 0.3978 * len2;

  return len;
}

２次元のベクトル(x,y)の長さを求めます。
「sqrt(x*x+y*y)」の近似となります。
これも、速度としてはまあまあ速いのですが、「誤差が大きい」ところに注意してください。
この関数で、たとえば円を描くと円周がカクカクします。
でも、アタリを求めるのには使えるでしょう。

※Windows環境では、なんと「sqrt(x*x+y*y)」のほうが上の関数よりも速いです。

mod：除算の余りを求める

double myMod(double x,double y)
{
  return (x - y * (double)((int)(x / y)));
}

xをyで割ったときの余りを求めます。yは0よりも大きい値を指定してください。
これは１行で計算できますので、インライン化したりすると、「fmod」よりも速いです。
特に、OS Xでは劇的な差がありますので、余りを求めたい場合はこの方法で行くと良いです。

実数の小数部だけを求めたい、という場合は「x - (double)((int)(x / y))」とするといいですね。

以上で、使えるとすると逆数平方根を求める関数でしょうか。これがあれば、ベクトルの単位ベクトルを求めるのなどは高速化につなげることができます。
余りを求める関数も（あまり使いどころはないかもしれないですが）こちらのほうがいいでしょう。
また、他にナイスなアルゴリズムをご存じの方はお教えいただければ幸いです。