平方根を使わないピタゴラス加算

二次元座標上の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()関数にも引けを取らない速度まで改善させることが出来ます。

また、オーバーフローのリスクが少ないこと、四則演算のみで確かに計算できていることにも注目してください。

実装の簡単さ・分かりやすさから、つい安直な方法を選んでしまいがちです。が、一時は面倒でも確実な方法を採ったほうが、最終的に面倒が少ないに違いないのです。

担当:田山(ピタゴラス加算の記号もあるらしいのですが、どう見ても排他的論理和です)