二次元座標上の2点間の距離を求めたいとき、複素数の絶対値を求めたいとき、その他いろいろなときに x = √(a2 + b2) という式を使います(これをpythagorean additionと言うそうです)。これを実装するとき、
sqrt( a*a + b*b ) |
という文がよく使われますが、この文を無闇に使っていると、あまり嬉しくない事態を引き起こすことになります。
さて、この文ののどこがいけないのでしょうか。
賢明な皆さまはもうお気づきでしょうが、a*a + b*bはオーバーフローを引き起こすおそれがあるのです。変数の値域を保証できるのならともかく、そうでなければ次に挙げる3つの対策のいずれかを採るべきです。
1.式変形
a' = max( |a|, |b| ),
b' = min( |a|, |b| )
とおくと、
√(a2 + b2) = a'√(1 + (b'/a')2)
と変形できます。この式を使うとオーバーフローのリスクを減らすことができます。
実装する際は、うっかり0除算のチェックを忘れないように注意してください。
2.hypot()関数を使う
C ( C++ )言語の math.h ( cmath ) ヘッダには hypot() という関数が定義されていて、double型の引数 a, b に対し √(a2 + b2) をdouble型で返してくれます。この関数は(返り値がdouble型に収まる限り)オーバーフローを起こさないように設計されているはずですので、安心して使うことができます。
3.Moler-Morrison Algorithmを使う
さて、こちらが今日の本題です。
Cleve Moler と Donald Morrison は、√(a2 + b2)を求めるための、短く、堅牢で、正確で、そして四則演算しか使わないアルゴリズムを発表しました。これをC++言語で書き直したものを御覧下さい。
double Molar_Morrison(double a, double b){ a = fabs(a); b = fabs(b); if (a < b) swap(a, b); if (b == 0) return a; const int ITERATION_NUMBER = 3; double s; for (int i=0; i<ITERATION_NUMBER; i++){ s=(b/a)*(b/a); s/=4+s; a+=2*a*s; b*=s; } return a; } |
このアルゴリズムの動作の肝を説明します。
√((a+2as)2 + (bs)2) は、式変形をしてみると √(a2 + b2) に全く等しいことが分かります。また明らかに s<1なので、bの値は0に収束していきます。したがって、aは次第に√(a2 + b2)に収束していくと言えるのです。
さて、プログラム中に ITERATION_NUMBER という定数が登場します。
これは「この回数だけループを繰り返せばdouble型で十分な精度が得られる」という数字です。ループを3回繰り返せば20桁の精度で解を求めることができ、また1回繰り返すごとに精度の桁数は少なくとも3倍になることが分かっています。
このプログラムは実行速度の面では若干の不安が残りますが、forループをアンローリングし、コンパイラの最適化オプションを使うなどの高速化を施せば、hypot()関数にも引けを取らない速度まで改善させることが出来ます。
また、オーバーフローのリスクが少ないこと、四則演算のみで確かに計算できていることにも注目してください。
実装の簡単さ・分かりやすさから、つい安直な方法を選んでしまいがちです。が、一時は面倒でも確実な方法を採ったほうが、最終的に面倒が少ないに違いないのです。
担当:田山(ピタゴラス加算の記号もあるらしいのですが、どう見ても排他的論理和です)