正の平方根を求めるアルゴリズム

はじめに

$\sqrt{a}$を求めます。基本的な考え方としては、関数$f\left(x\right) = x^2 - a $と$x$軸との正の共有点の$x$座標が$\sqrt{a}$であることを利用して、計算します。

１．平均を取って近似していく方法

関数 $f\left(x\right)$と$x$軸の交点を含む範囲の左端を$left$、右端を$right$とします。このとき、$left < \sqrt{a} < right$です。
$f\left(\frac{left + right}{2}\right) < 0$ のとき、$\frac{left + right}{2} < \sqrt{a} < right$ ・・・①
$f\left(\frac{left + right}{2}\right) > 0$ のとき、$left < \sqrt{a} < \frac{left + right}{2}$・・・②
①のとき、$left = \frac{left + right}{2}$ とすると、$left < \sqrt{a} < right$
②のとき、$right = \frac{left + right}{2}$ とすると、$left < \sqrt{a} < right$
この操作を繰り返すことで近似していきます。今回は結果の値として$right$を返します。

sqrt1.cpp

double sqrt1(double a) {
    double left = 0.0;
    double right = a;
    double mid = 0.0;
    for(int i = 0; i < 1000; i++) {
        mid = (left + right) / 2;
        if(mid * mid - a < 0) left = mid;
        else right = mid;
    }
    return right;
}

2. 微分を用いて近似していく方法

$\sqrt{a} < b$ を満たす実数$b$を考えたとき、点$(b,\, f\left(b\right))$における関数$f\left(x\right)$の接線を$l$とします。この接線$l$は、$f'\left(x\right) = 2x$より、
$l: y = 2b(x - b) + b^2 - a = 2bx - (b^2 + a)$
と表されます。ここで、接線$l$と$x$軸との交点の$x$座標は、$2bx - (b^2 + a) = 0$より、$x = \frac{b^2 + a}{2b}$となります。ここで、$\lim_{b \to \sqrt{a}}{x} = \frac{2a}{2\sqrt{a}} = \sqrt{a}$であり、また、$x < b$であることから、$b = x = \frac{b^2 + a}{2b}$として操作を繰り返していけば、$b$が$\sqrt{a}$に近づくことがわかります。

sqrt2.cpp

double sqrt2(double a) {
    double b = a;
    for(int i = 0; i < 1000; i++) {
        b = (b * b + a) / (2 * b);
    }
    return b;
}

おまけ

計算量について

この2つの方法は、どちらの方が先に$\sqrt{a}$に近づくのか、求めたい精度の有効桁数を$N$として求めてみましょう。
まず、1番目の方法は、$a × 10^{N-1}$個の数の中から1つを選ぶ二分探索と同じですので、計算量は$O(\log(a×10^{N-1})) → O(N)$です。
それに対して、2番目の方法は簡単には分かりません。そこで、以下のコードを実行しました。

sqrt3.cpp

#include <cstdio>
#include <unistd.h>

double sqrt3(double a) {
    double b = a;
    for(int i = 0; i < 10; i++) {
        b = (b * b + a) / (2 * b);
        printf("step:%d val:%.20f\n", i + 1, b);
        usleep(100000);
    }
    return b;
}

int main() {
    sprt3(2.0);
    return 0;
}

すると、以下のような実行結果が得られました。

step:1 val:1.50000000000000000000
step:2 val:1.41666666666666674068
step:3 val:1.41421568627450988664
step:4 val:1.41421356237468986983
step:5 val:1.41421356237309492343
step:6 val:1.41421356237309492343
step:7 val:1.41421356237309492343
step:8 val:1.41421356237309492343
step:9 val:1.41421356237309492343
step:10 val:1.41421356237309492343

この結果をまとめると、
$N=1$のとき、$step=1$
$N=2,3$のとき、$step=2$
$N=4,5,6$のとき、$step=3$
$N=7〜12$のとき、$step=4$
…と続いていくことから、$O(\log N)$であると推測できます。

もちろん1番目の場合は$a$の値(というより$right$を初期化する際の値)によって、2番目の場合は$b$を初期化する際の値によって計算に必要な時間は変化します。そのため一概には言えませんが、一般的に2番目の方法の方が優れたアルゴリズムであると考えられるでしょう。

追記

2番目の方法がニュートン法であることに後から気付きました。とても有名なアルゴリズムです。何で気付かなかったんだろ。ちなみに1番目の方法も二分法ってやつです。これまた有名なアルゴリズムです。