2010-11-20

平方根のアルゴリズム

平方根である。sqrtである。私の数学知識は非常に乏しく、かろうじて平方根の何たるかを解するばかりであるが、ふと、平方根を計算するアルゴリズムが知りたくなった。理由は、constexprなsqrtを実装してみたくなったからだ。

日本語では、いい情報がWeb上に存在しないので、ここに書いておく次第。この説明は、私のように数学の知識を持たない人間にも理解できるはずである。

ここで使うのは、バビロニア人の方法(Babylonian method)と呼ばれている、歴史あるアルゴリズムである。この方法は、手動でも、プログラミングでも、非常に簡単に計算できる、大変便利で汎用的なアルゴリズムである。しかも、速度もそれほど遅くない。

√Sに対して、

  1. 任意の正の整数を初期値X0と定める(できるだけ平方根に近い値が望ましい)
  2. Xn+1を、Xnと S / Xnの平均値とする(平均値は相加平均で求める)
  3. 必要な精度が得られるまで、ステップ2を繰り返す

では、さっそく手動で計算してみよう。今、√123を求めたいとする。

上のアルゴリズムに当てはめると、S = 123である。x0は、とりあえず1とする。すると

x1 = ( x0 + 123 / x0 ) / 2 = ( 1 + 123 / 1 ) / 2 = 62
x2 = ( 62 + 123 / 62 ) / 2 = 31.9919
x3 = 17.9183
x4 = 12.3914
x5 = 11.1588
x6 = 11.0905

たったの6回の計算で、それなりの精度の平方根を求められる。

さて、これをコードに落としこむ際に、初期値をどう決めるかという問題がある。アルゴリズム的には、正の整数でさえあれば、正しい答えは出る。しかし、ここで条件分岐などを用いて複雑な方法で初期値を決めてしまっては、本末転倒である。大抵のsqrtの実装では、できるだけ簡単な値が選ばれる。たとえば、Sの二分の一など。次に、何度ステップを繰り返すかということであるが、これは、そもそも浮動小数点数の性質からして、値を正確に表すことができないので、単に、前回の値を比較するという方法を取る。つまり、値を比較して差がわからなければ、もうそれ以上ステップを踏む必要がないということである。

namespace hito {

double sqrt( double s )
{
    double x = s / 2.0 ; // Is there any better way to determine initial value?
    double last_x = 0.0 ; // the value one before the last step

    while ( x != last_x ) // until the difference is not significant
    { // apply Babylonian method step
        last_x = x ;
        x = (x + s / x) / 2.0 ;
    }
    
    return x ;
}

}

さて、この関数をconstexpr関数にするのは、簡単である。単に副作用をなくし、ループを再帰で置き換えればいいのだ。

namespace hito
{
    constexpr double sqrt_impl( double s, double x, double last )
    {
        return x != last ? sqrt_impl( s, (x + s / x) / 2.0, x ) : x ;
    }

    constexpr double sqrt( double s )
    {
        return sqrt_impl( s, s / 2.0, 0.0 ) ;
    }
}

ただし、コンパイル時に平方根を計算したい需要がそれほどあるとも思われない。

さらに詳しく知りたい人のために、有益なサイトを挙げておく。

根本的な計算方法は、Wikipediaに詳しく載っている。

Methods of computing square roots - Wikipedia, the free encyclopedia

具体的なコードで読みたいという人のためには、C++によるコード集もある。ちなみに、この記事は素晴らしい出来で、各アルゴリズムの出典まで明記している。

Best Square Root Method - Algorithm - Function (Precision VS Speed) - CodeProject

予想していた通り、浮動小数点数の内部表現をあてにしたビット演算のトリックがかなりある。中でも興味深いのは、謎のマジックナンバーで減算した後、手動のニュートン法を二回分行うコードである。噂によると、Quake 3を開発した伝説のゲームプログラマー、John Carmackによって発明されたのだと言われている。

6 comments:

Akira said...

その収束条件だとまれに収束しない場合があります。浮動小数点数で表現可能な値のうち√Sの真値に最も近い前後の値で振動し続けます。
その場合、constexprだとコンパイル時にstack overflowが発生するのでしょうか?

Akira said...

と思ったんですが、収束しない場合があるのは一般的なニュートン法の話で、

x = (x + s / x) / 2.0 ;

この計算なら大丈夫そうですね。

江添亮 said...

そう、大丈夫なんですよ。すごいですよね。
Babylonian methodの文献上の初出は、アレクサンドリアのヘロンだそうです。約二千年前です。

それがいまだに使えるんだから。

江添亮 said...

ただ、ニュートン法の派生物と考えることもできるそうで。
ただ、歴史的には、ニュートン法の考案はだいぶ最近の話ですが。

Anonymous said...

なかなか興味深いアルゴリズムですね.
実数では収束するが四元数などでは変な初期値から始めると収束しないだろうなあと思い,実際にプログラミングして見た.変な初期値ではだめ.
挟み撃ちアルゴリズムの性格があるから大き目の初期値から行くと収束する.
その反面収束が遅くなる.例:Q=(10,2,2,2)では,
Q^0.5=(3.20804,0.311717,0.311717,0.311717)になったのは50回もかかった.

Anonymous said...

可除環の最後8元数でもこのアルゴリズムは可能です.実際にコーディングしてみました.
収束性は4元数と同じレベルです.Babylonian methodの本質は除算可能性です.