平方根の値を計算する方法(二分法・ニュートン法・開平法)について

 \sqrt{5} = 2.2360679 \dots であることが知られていて「富士山麓オウム鳴く」のごろ合わせも有名です。実際に 2.2360679^{2} = 499999965341041 \lt 5 \lt 500000010062400 = 2.2360680^{2} であることから小数第6,7位が80ではなく79が正しいことが確認できます。

ところで平方根の値を計算する効率の良い方法を知っていますか?

この記事では平方根の値を求める方法として、計算方法を3つ紹介し、PythonとRustでの実装を示します。

  1. 二分法
  2. ニュートン法
  3. 開平法

二分法

求めたい平方根の値は  0 \le \sqrt{x} \le max(1.0, x) の範囲にあるため、二分法で範囲を絞っていくことで近似値を求めることができます。

1回のループごとに解の存在範囲が 1/2 倍になります。浮動小数点数型を使う場合は型の精度に気を付けましょう。

def sqrt_by_bisection(k, epslilon=0.001):
    """k-epsilon <= x*x < k であるxを返す。"""
    upper = 1 if k < 1 else k
    bottom = 0
    while upper - bottom >= epslilon:
        m = bottom + (upper - bottom) / 2
        if m**2 > k:
            upper = m-epslilon
        else:
            bottom = m
    return bottom
// k-epsilon <= x*x < k であるxを返す。
fn sqrt_by_bisection(k: f64, epsilon: f64) -> f64 {
    debug_assert!(k >= 0.0);
    let mut upper = if k < 1.0 { 1.0 } else { k };
    let mut bottom = 0.0;
    while upper - bottom >= epsilon {
        let m = bottom + (upper - bottom) / 2.0;
        if m * m > k {
            upper = m - epsilon;
        } else {
            bottom = m;
        }
    }
    bottom
}

ニュートン法

 \sqrt{x} の値を求めるという問題は、 f(x) = x^{2} - k として f(x) = 0 となる  x を求める問題と言い換えることができ、この問題はニュートン法で解の近似値を求めることができます。

この方法は二分法よりも収束が速いことが知られています。

2次収束と呼ばれ,実際に計算してみるとわかるように,非常に速い. 大ざっぱにいって「誤差が毎回2乗で減る」あるいは「正しい桁数が反復1回ごとに倍に増える」といえる.

『数値計算の常識』P60より

def sqrt_by_newton(k, epsilon=0.001):
    def f(x): return x - (x * x - k) / (2.0 * x)
    if k == 0.0:
        return 0.0
    x = k
    next_x = f(x)
    while abs(x - next_x) >= epsilon:
        x = next_x
        next_x = f(x)
    return next_x
fn sqrt_by_newton(k: f64, epsilon: f64) -> f64 {
    let f = |x| x - (x * x - k) / (2.0 * x);
    debug_assert!(k >= 0.0);
    if k == 0.0 {
        return 0.0;
    }

    let mut x = k;
    let mut next_x = f(x);
    while (x - next_x).abs() >= epsilon {
        x = next_x;
        next_x = f(x);
    }

    next_x
}

実測:ニュートン法は二分法よりも収束が速いのか

ニュートン法のほうが速いし実装も簡単です。

上記のPython実装を使って精度を変更しながらループ回数との関係を調べました。浮動小数点数型の誤差の問題でi=14以降計算することが出来ませんでした。Decimalモジュールを使うとより多い桁数を正しく求めることができるかも知れません。

f:id:matsu7874:20210109035543p:plain

開平法

手計算(筆算)で行いやすい方法として開平法という方法があります。

 10a + b \le \sqrt{x} であるような  b (10a + b)^{2} = 100a^{2} + (2 \times 10a + b)bという形に変形して求めていく方法です。詳しい方法は参考のリンクを参照してください。

実装上、整数で演算するので多倍長整数が使える言語(Pythonなど)ではかんたんに多くの桁数を正確に計算することが出来ます。

def sqrt_by_extraction(n, prec=16):
    chunks = []
    k = n
    while k > 0:
        chunks.append(k % 100)
        k //= 100
    sqrt = []
    a = 0
    b = 0
    for _ in range(prec):
        if chunks:
            c = chunks.pop()
        else:
            c = 0
        a *= 10
        b = b * 100 + c
        for i in range(10)[::-1]:
            t = (a + i) * i
            if t <= b:
                b -= t
                a += 2 * i
                sqrt.append(i)
                break
    return sqrt
fn sqrt_by_extraction(n: i64, prec: usize) -> Vec<i64> {
    debug_assert!(n >= 0);
    debug_assert!(prec <= 16);
    let mut chunks = vec![];
    let mut k = n;
    while k > 0 {
        chunks.push(k % 100);
        k /= 100;
    }
    eprintln!("{:?}", chunks);
    let mut sqrt = vec![];
    let mut a = 0;
    let mut b = 0;
    for _ in 0..prec {
        let c = chunks.pop().unwrap_or(0);
        a *= 10;
        b = b * 100 + c;
        for i in (0..10).rev() {
            let t = (a + i) * i;
            if t <= b {
                b -= t;
                a += 2 * i;
                sqrt.push(i as i64);
                break;
            }
        }
    }

    sqrt
}

まとめ

平方根の値を計算する3種類の方法(二分法・ニュートン法・開平法)について紹介し、それぞれPythonとRustでの実装例を示しました。 大きな桁数が必要であれば開平法を、そうでなければ高速なニュートン法を使うと良さそうです。