平方根のアルゴリズム

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

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

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

√Sに対して、

1. 任意の正の整数を初期値X₀と定める（できるだけ平方根に近い値が望ましい）
2. X_n+1を、X_nと S / X_nの平均値とする（平均値は相加平均で求める）
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によって発明されたのだと言われている。