1. はじめに
バビロニアの平方根は、ニュートン法の特殊系として知られており、平方根を効率よく近似出来るアルゴリズムである。
この記事では、平方根を高速に求める手法であるバビロニアの平方根について説明し、Pythonで実装を行う。また、二分探索での実装と、組み込み関数である pow とのベンチマークを比較し、その性能を検証する。
2. ニュートン法
バビロニアの平方根について説明する前に、ニュートン法について解説をする。これは、バビロニアの平方根がニュートン法の特殊系であるためである。
ニュートン法は、方程式の解を高速に計算するアルゴリズムとして知られている。以下の手順で、解の近似を進めていく。
- 方程式を $f(x) = 0$ に変形する
- 適当な初期値 $x_0$ を与える
- 以下の反復式に従って、$x_{n+1}$ を計算する： $$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$
- $x_{n+1}$ が十分に収束するまで、3 を繰り返す
この手法は初期値が解に十分近い場合、非常に高速に収束するという特徴を持つ。
ここで、$\sqrt{2}$を例にとりニュートン法を実践していく。
- 方程式を $f(x) = 0$ に変形する 解を$x$とする。$x = \sqrt{2}$より、両辺を2乗して$x ^ 2 - 2 = 0$ として$f(x) = 0$の形に帰着できた。
- 適当な初期値 $x_0$ を与える 次に初期値 $x_0$ を取る。一般的に、平方根の初期値としては対象の値の半分程度が選ばれることが多い。ここでは$x_0 = 1.0$を初期値として用いる。
- 反復式 $$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$ を使って$x_{n+1}$を計算していく。これは$x_{n+1}$が十分に収束するまで繰り返される。 $f(x_n) = x_n ^ 2 - 2$、$f'(x_n) = 2 x_n$、$x_0 = 1.0$ として計算をしていく。これを計算していくと以下のようになる。
- グループA：二分探索とバビロニアの平方根
- グループB：上記2手法に加えて、組み込み関数である
powを用いた平方根計算
1.0 # 1回目の試行
1.50000000 # 2回目の試行
1.41666667 # 3回目の試行
1.41421569 # 4回目の試行
1.41421356 # 5回目の試行
桁数は小数点第8桁まで求めた。この時5回の試行で、十分な精度の近似をすることができた。
3. バビロニアの平方根
ここでは、題名にもある通りバビロニアの平方根について説明をする。前段でも述べた通り、バビロニアの平方根はニュートン法の特殊系である。
求めたい平方根の値を$S$と置く。これにニュートン法を適用するには、次のように関数を定義する。
$$
f(x) = x ^ 2 - S
$$
この時の導関数は、
$$
f'(x) = 2x
$$
これらをニュートン法の反復式に代入すると、
$$
x_{n+1} = x_n - \frac{x_n^2 - S}{2x_n}
$$
これを整理して、
\begin{align*}
x_{n+1} &= x_n - \frac{x_n^2 - S}{2x_n} \\
&= x_n - \frac{x_n}{2} + \frac{S}{2x_n} \\
&= \frac{x_n}{2} + \frac{S}{2x_n} \\
&= \frac{1}{2} (x_n + \frac{S}{x_n})
\end{align*}
この式が、いわゆるバビロニアの平方根である。ニュートン法の特殊系であるため、わざわざこちらを使わずとも、ニュートン法を使えば同じ程度の速度で平方根を求めることは可能である。しかし、こちらの式の方が単純であり、浮動小数点演算との相性も良いという利点がある。
4. ベンチマークの計測
ここまでは、ニュートン法と、その派生形であるバビロニアの平方根について説明をしてきた。ここでは、ニュートン法の実行速度を他の手法と比較して、性能を検証する。
実行速度の計測には、平方根の計算を 10000 回実行し、1度の計算ごとに、最も高速であった手法に +1、最も低速であった手法に -1 の得点を与える。すべての試行を通じて最も得点が高かった手法を、最も高速であると判断する。
平方根の計算対象には、1から100000までの整数の中からランダムに選ばれた値を用いる。各手法においては、浮動小数点計算の誤差を考慮し、十分な収束とみなすための許容誤差 $\alpha = 10^{-9}$ を設定している。反復は、この誤差を下回るまで繰り返される。
実行速度の比較は、次の2つのグループに分けて行う。
4.1. グループAの実行結果
グループAの実行結果を、下記に示す。
| 手法 | 勝ち数 | 負け数 | 勝率 |
|---|---|---|---|
| 二分探索 | 258 | 9731 | 2.58 % |
| バビロニアの平方根 | 9731 | -258 | 97.31 % |
| 引き分け | 3 | - | - |
実行結果から、今回の計測の条件においては必ずしもバビロニアの平方根が常に最速であるとは限らないことが分かる。しかし、全体の勝率は97.31 %に達しており、二分探索と比較して圧倒的に高速であることが分かる。
4.2. グループBの実行結果
グループBの実行結果を、下記に示す。
| 手法 | 勝ち数 | 負け数 | 勝率 |
|---|---|---|---|
| 二分探索 | 5 | -9650 | 0.05 % |
| バビロニアの平方根 | 112 | -300 | 1.12 % |
組み込み関数 pow
|
9883 | -50 | 98.83 % |
| 引き分け | 0 | - | - |
実行結果から、二分探索とバビロニアの平方根と比較して、組み込み関数 pow が圧倒的に高速であることがわかる。全体の勝率は98.83 %に達しており、通常使用の場合は pow を使うことが推奨される。
5. Python コード
前段のベンチマークの計測で使用したPythonコードを下記に示す。
import time
import random
"""
実行時間の計測の為に使用する `measure_time` は、デコレータを用いて実装している。
"""
def measure_time(func):
def timer(*args):
start = time.time()
result = func(*args)
end = time.time()
total_time = end - start
print(f"関数 {func.__name__} の実行時間: {total_time:.6f} 秒")
return result, total_time
return timer
@measure_time
def sqrt_babylonian(n):
guess = n / 2
alpha = pow(10, -9)
while abs(guess * guess - n) > alpha:
guess = (guess + n / guess) / 2
return f"{guess:.6f}"
@measure_time
def sqrt_binary(n):
alpha = pow(10, -9)
left = 0
right = n
while right - left > alpha:
mid = left + (right - left) / 2
if mid * mid < n:
left = mid
else:
right = mid
return f"{left:.6f}"
@measure_time
def sqrt_pow(n):
return pow(n, 0.5)
def compare_speeds(name1, time1, name2, time2):
print("実行時間の差: ", end="")
if time1 < time2:
print(f"{name1} の方が {time2 - time1:.6f} 秒早い")
else:
print(f"{name2} の方が {time1 - time2:.6f} 秒早い")
def benchmark_once(num):
res_bin, time_bin = sqrt_binary(num)
res_new, time_new = sqrt_babylonian(num)
res_pow, time_pow = sqrt_pow(num)
print(res_bin)
print(res_new)
print(res_pow)
compare_speeds("sqrt_babylonian", time_new, "sqrt_binary", time_bin)
compare_speeds("sqrt_babylonian", time_new, "sqrt_pow", time_pow)
compare_speeds("sqrt_binary", time_bin, "sqrt_pow", time_pow)
timings = {
"sqrt_binary": time_bin,
"sqrt_babylonian": time_new,
"sqrt_pow": time_pow,
}
fastest = min(timings, key=timings.get)
slowest = max(timings, key=timings.get)
return fastest, slowest
def main():
COUNT = 10000 # 試行回数
# [a, b] => a: 勝利数, b: 敗北数
rank = {
"sqrt_binary": [0, 0],
"sqrt_babylonian": [0, 0],
"sqrt_pow": [0, 0]
}
ties = 0 # 引き分けの回数
for _ in range(COUNT):
num = random.randint(1, 100000) # 求める平方根を乱数生成によって取得する。
fastest, slowest = benchmark_once(num)
if fastest == slowest:
ties += 1
continue
rank[fastest][0] += 1
rank[slowest][1] -= 1
print(f"COUNT: {COUNT} 回")
print(f"ties: {ties} 回（同率試行）")
print(rank)
if __name__ == "__main__":
main()
6. 最後に
今回は、平方根を高速に求める手法であるバビロニアの平方根について解説し、その実装とベンチマークの計測を行った。
ベンチマークの結果から、バビロニアの平方根は二分探索よりも高速であり、平方根を求めるアルゴリズムとして有力な選択肢であることが分かった。
一方で、実用上は組み込み関数（pow など）を使用する方が、圧倒的に高速であることも明らかとなった。アルゴリズムの学習や理解を目的とする場合を除けば、基本的には組み込み関数の使用が推奨される。
7. 参考資料
[1]. 高校数学の美しい物語 - ニュートン法の解説とそれを背景とする入試問題