Pythonで平方根を計算する

Pythonを使って平方根を計算してみます。といっても、x**(0.5)とかsqrt関数で計算しようという話ではありません。頭の体操で、手計算に近い方法で計算しようという趣旨です。

ごくごく単純な方法

中学校程度の数学の知識で計算する、とてもシンプルな方法です。

#1 原始的な方法で平方根を求める
def sqrt_1(num,prec):
    x=0
    while True:
        if num < (x+1)**2:
            break
        x+=1    
    for digit in range(1,prec+1):
        i=0
        while True:
            if num <(x+10**(-digit)*(i+1))**2:
                break
            i+=1
        x+=10**(-digit)*i
    return round(x,prec)

sqrt_1(2,8)

1.41421356

関数sqrt_1に数値（num:小数でも可）と求める桁数（prec）を渡します。4行目のwhileループではxを、その2乗がnumより大きくなるまで0から1ずつ増やしていきます。この結果、xの値はnumの平方根の整数部分になります。

次に小数部分を7行目のforループの中で、1からprecまで、1桁ずつ計算していきます。ここでは、小数点第1位を例にするとxに、やはりその2乗がnumより大きくなるまで0.1ずつ増やしていきます。ここで変数degitが小数点以下の桁数を示していて、10の-digit乗ずつ増やすことにより、求める桁数を増やすことができます。このことから最後の桁での値は切り捨てた値となることに注意してください。非常にシンプルな方法ですが、ルート2は1.41421356と正しく計算されました。

ニュートン法による方法

ニュートン法は、初期値としてとりあえず答えより少し大きな値を初期値$x_0$として与えておき、そこから微分の考え方を使い徐々に真の値に近づけていく方法です。そこで、次の式で計算します。

$x_{k+1}=x_{k}-\displaystyle\frac{f(x_k)}{f^{\prime}(x_k)}$

numの平方根を求める場合、$f(x)=x^2-num$となるので、

$x_{x+1}=x_{k}-\displaystyle\frac{{x_{k}}^2-num}{2 x_k}=\frac{1}{2}\left(x_k+\frac{num}{x_k}\right)$

#2 ニュートン法 
num=2
x=1.5
for i in range(5):
    x=(1/2)*(x+num/x)
    print(f'{x:.20f}')

1.41666666666666651864
1.41421568627450966460
1.41421356237468986983
1.41421356237309492343
1.41421356237309492343

初期値を1.5として小数点20桁までを計算しています。正しい結果は1.4142135623730950488016887なので、4回ループさせるだけで小数点第14位まで正しく計算されます。

decimalモジュールを使いもう少し正確に計算する

上記の方法で小数点15桁目以降で差異が生じていますが、Pythonで小数点以下の計算をすると誤差が生じることからニュートン法による問題ではないように思われます。そこでdecimalモジュールを使い100桁まで計算してみます。

#3 ニュートン法 
from decimal import * 
getcontext().prec = 100
num=2
x=1.5
for i in range(8):
    x=(Decimal('1')/Decimal('2'))*(Decimal(str(x))+Decimal(str(num))/Decimal(str(x)))
    print(f'{x:.100f}')

1.4166666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666660
1.4142156862745098039215686274509803921568627450980392156862745098039215686274509803921568627450980390
1.4142135623746899106262955788901349101165596221157440445849050192000543718353892683589900431576443400
1.4142135623730950488016896235025302436149819257761974284982894986231958242289236217849418367358303560
1.4142135623730950488016887242096980785696718753772340015610131331132652556303399785317871612507104750
1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276416020
1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415720
1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415720

正しくは1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573なので100桁程度は正しく計算できることになります。たかだか8回程度でこの制度ですから恐るべしといわざるを得ません。

ニュートン法で求める桁数を制御する

ニュートン法では少ないループで平方根をかなり良い精度で計算することができますが、ここでは小数点以下20位までを正確に計算する方法を考えます。よく取られる方法として$x_{k+1}$と$x_{k}$を比較して、その差が$10^-20^に収まっていればループを終わらせるという方法があります。

from decimal import *
e=20
getcontext().prec = e
num=2
x=1.5

cnt=0
while True:
    cnt+=1
    x1=(Decimal('1')/Decimal('2'))*(Decimal(str(x))+Decimal(str(num))/Decimal(str(x)))
    if abs(Decimal(x1)-Decimal(x))<Decimal('0.1')**Decimal(str(e)):
        break
    x=x1    
print(f'{e}桁の精度：{cnt}回のループで{x1}')

20桁の精度：5回のループで1.4142135623730950488

精度の計算をするとdecimalモジュールを使うなど、デリケートなものです。