さて、この文ののどこがいけないのでしょうか。
賢明な皆さまはもうお気づきでしょうが、a*a + b*bはオーバーフローを引き起こすおそれがあるのです。変数の値域を保証できるのならともかく、そうでなければ次に挙げる3つの対策のいずれかを採るべきです。

1.式変形

a' = max( |a|, |b| ),
b' = min( |a|, |b| )
とおくと、

√(a² + b²) = a'√(1 + (b'/a')²)

と変形できます。この式を使うとオーバーフローのリスクを減らすことができます。
実装する際は、うっかり0除算のチェックを忘れないように注意してください。

2.hypot()関数を使う

C ( C++ )言語の math.h ( cmath ) ヘッダには hypot() という関数が定義されていて、double型の引数 a, b に対し √(a² + b²) をdouble型で返してくれます。この関数は（返り値がdouble型に収まる限り）オーバーフローを起こさないように設計されているはずですので、安心して使うことができます。

3.Moler-Morrison Algorithmを使う

さて、こちらが今日の本題です。
Cleve Moler と Donald Morrison は、√(a² + b²)を求めるための、短く、堅牢で、正確で、そして四則演算しか使わないアルゴリズムを発表しました。これを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)² + (bs)²) は、式変形をしてみると √(a² + b²) に全く等しいことが分かります。また明らかに s<1なので、bの値は0に収束していきます。したがって、aは次第に√(a² + b²)に収束していくと言えるのです。

さて、プログラム中に ITERATION_NUMBER という定数が登場します。
これは「この回数だけループを繰り返せばdouble型で十分な精度が得られる」という数字です。ループを3回繰り返せば20桁の精度で解を求めることができ、また1回繰り返すごとに精度の桁数は少なくとも3倍になることが分かっています。

このプログラムは実行速度の面では若干の不安が残りますが、forループをアンローリングし、コンパイラの最適化オプションを使うなどの高速化を施せば、hypot()関数にも引けを取らない速度まで改善させることが出来ます。

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

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

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