平方根を使わずに高速で２点間の距離を近似する

２点間の距離の計算では平方根が必要になりますが、平方根は少し重い計算です。ということで、平方根を使わず、掛け算・割り算・足し算と絶対値・最大・最小だけで距離を近似する方法についての記事を翻訳してみました。
flipcode - Fast Approximate Distance Functions
(12:02 補足：おそらく今の標準的なCPUでやる意味はほとんどないと思います。近似のアプローチとして面白いというくらいの話。Z80でやりましょう)

距離関数高速近似

by Rafael Baptista (27 June 2003)

２点間のユークリッド距離を求める計算式は次のようになる。

二次元では次のようになる。

この関数の計算には、平方根が必要になる。これは最近のコンピュータでも高価な計算である。平方根は逐次近似によって求められる。つまり、コンピュータは平方根近似のループを行って、与えられた許容誤差に収まるまで繰り返す。もし、もっと制約の多いハードウェアで作業していて、ハードウェアでの平方根実装がなければ、小さい数の平方根の計算でさえも禁止されるだろう。

多くの計算で、二乗距離での比較が可能で、平方根をまったく使わないようにできる。例えば、必要な計算が距離の比較だけだと考えられるなら、二乗距離での計算は等価である。しかしながら、本当に距離を求めないといけない計算もある。あげてみると、ベクトルの正規化だとか、球面のバウンディングボリュームを使った衝突システムの実装だとか。

これらの場合に、標準的な距離関数を計算する余裕がないとして、距離関数のちょっといい近似をして、計算の簡単な線形演算だけで構成される関数クラスがある。理論上はもっと演算を追加することで、線形近似を保ったまま任意の精度で非線形関数を計算できる。実際これは、コンピュータで実装されるときに平方根関数が行うことと似ていないわけではない。ここで述べようとしている関数は、概して距離を2.5%の平均誤差で近似し、少ない線形計算の場合でだいたい最大5%の誤差あるが、そのため実行が非常に高速である。最大誤差と平均誤差のトレードオフの調整は可能で、そして誤差を適切な値にもってくるように関数を構成することも可能である。

最初に、正しい距離関数の3Dグラフを見てみよう。この3Dグラフはxとyで-1から1までの範囲のユークリッド距離をあらわす。グラフの高さは負の距離を表す。xとyが0のとき、距離は0であり、円錐の頂上になる。xとyの値を増やすと、それに応じて円錐の点は下がっていく。このグラフの断面は円のようになる。

ここに、xとyの絶対値の最大と最小の、xとyの範囲が-1から1までのグラフがある。最大のグラフは4面のピラミッドになる。xとyが0になる点は、このピラミッドの頂点になる。この断面は正方形になる。

最小値のグラフはピラミッドとは逆の操作になることがわかる。その断面は十字になる。

その結果、このふたつの関数の線形合成を行うと、距離関数へのちょっとよい近似を得ることができる。 最大と最小の異なった比率の値は異なる属性の近似を生成する。ここに挙げたものは、最大誤差を抑えるように調整してある。この関数を平均誤差が最小になるよう、もしくは平均二乗誤差が最小になるように調整することもできる。最大誤差になるのはxとyが0に近づくところであることがわかる。そして、これらは断面の形が理想円からもっともはずれる場所である。 (比率の項を調整することで、この点での誤差をゼロに減らすことは可能である。しかし、ほかの場所での精度を大きく犠牲にする)

これらの点での誤差を補正するため特別に設計した項を追加することで、この関数の精度をあげることができる。この場合、xやyがより大きいところ、xとyの小さいほうの16倍になるところだけで関数の比率を変えればよい。この新しい関数のグラフには、xやyが0に近づくところにギザギザがややできあがっていることを見ることができる。目的にあった誤差になるまで、処理を繰り返して総誤差を減らすことができる。もちろん、線形の項を追加すると、関数の計算はより高価になる。

計算したい、xとyの絶対値の最大と最小を伴った距離関数の近似を証明することは簡単だ。正しい距離関数を考えると、与えられたどんな距離Xについても半径Xの円を得る。本質的には円のグラフを直線で近似することを試す。xとyの絶対値を取るのは、距離関数の結果を変えないことに注意する。最初にxとyの絶対値を最初に取ることで、問題を円のひとつの象限だけの近似に減らすことができる。もうすこし考えると、xとyの最大最小をとることで、問題を８分円だけに限定できる。これは、ちょっとよい近似を少数の直線だけで得るのに十分小さい弧である。

これらの関数は非常に簡単に計算でき、最近のハードウェアでは定数時間で計算することができる。今回使った係数は1024の分数として表せることに注意する。これはこの関数を整数レジスタと少しの掛け算だけを使って実装するできることを意味し、最終的にシフトを使うだけで結果をスケールダウンできる。2の乗数になる分母を使った係数を選べば、これが可能になる。どれだけ多くのビットを使うかは、レジスタサイズと必要な精度に依存する。

これはひとつの近似関数の実装例である。

u32 approx_distance( s32 dx, s32 dy )
{
   u32 min, max, approx;

   if ( dx < 0 ) dx = -dx;
   if ( dy < 0 ) dy = -dy;

   if ( dx < dy )
   {
      min = dx;
      max = dy;
   } else {
      min = dy;
      max = dx;
   }

   approx = ( max * 1007 ) + ( min * 441 );
   if ( max < ( min << 4 ))
      approx -= ( max * 40 );

   // add 512 for proper rounding
   return (( approx + 512 ) >> 10 );
} 

この近似は、ハードウェアが制限されているのであれば、掛け算さえ使わずに実装することも可能である。

u32 approx_distance( s32 dx, s32 dy )
{
   u32 min, max;

   if ( dx < 0 ) dx = -dx;
   if ( dy < 0 ) dy = -dy;

   if ( dx < dy )
   {
      min = dx;
      max = dy;
   } else {
      min = dy;
      max = dx;
   }

   // coefficients equivalent to ( 123/128 * max ) and ( 51/128 * min )
   return ((( max << 8 ) + ( max << 3 ) - ( max << 4 ) - ( max << 1 ) +
            ( min << 7 ) - ( min << 5 ) + ( min << 3 ) - ( min << 1 )) >> 8 );
} 

(翻訳ここまで)

単語

原文で出てきた単語について

cross section	断面
ideal	理想の
diviate	基準からはずれる
quadrant	象限
demonstrate	証明
constrain	制限
coefficient	係数
denominator	分母

参考になる本

こういった、実数値を求めるような計算を数値計算だとか数値解析というけど、機械学習とかやろうと思うと知っておく必要がある。
数値計算の入門としてはこの本がおすすめ。
扱う範囲がしぼってあって、読みやすい。例はBASICだけど、なんとなくわかると思う。

数値計算のはなし―パソコンを使って解こう

具体的な計算の方法は、定番のニューメリカルレシピ。まあ古くていろいろ批判があってバグもあって、本当は新しい版の訳がほしいところだけど、手始めに。本気でやるなら原著の最新版のほうがいいです。

ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ

あと、テイラー展開とかは知っておく必要があるので、無印数学ガールなど読んでおくといいと思います。

数学ガール (数学ガールシリーズ 1)