自作関数シリーズ第1回 ルート関数 (C言語、Java、Python、Ruby)

このシリーズは普段ライブラリを利用して求めているような計算を自作するとしたらどのようなコードになるのかを紹介するシリーズです。といってもOpenCVやBeautiful Soupのような大規模なライブラリの関数を自作するわけではありません。大規模なライブラリは眺めてこんな感じ何だなって理解するに留めて、素直に作ってくれた人に感謝して利用しましょう。さて記念すべき第1回はルートを求める関数を自作したいと思います。ちなみに普通にルートを求めるには

C言語: include<math.h>をしてからsqrt関数を使う

Java: Mathクラスのsqrtメソッドを使う

Python: import mathをしてsqrt関数を使う

Ruby: Math.sqrtを使う

今回ルートを求めるのに使う方法はニュートン法です。

ニュートン法は「f(x) = 0になるような値Xを探す時、ある値Xnにおける接線の切片Xn+1は、元の値Xnより真の値Xに近くなる。」という考え方に基づいています。

この考え方を数式に直すとXn+1とXnの漸化式は

Xn+1 = Xn - f(Xn)/f'(Xn)

となります。

下の図を例に見てみると確かにXn+1はXnよりも真の値Xに近いことがわかります。

ニュートン法の数学的な説明をもっと詳しく知りたい人は

wikipediaなどを参考に読んでみてください。

https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E6%B3%95

さて今回の本題ルートを上で紹介したニュートン法で求めていきましょう。

x = √a

と置き式変形すると

x^2 - a = 0

ここで左辺をf(x)と置くと

f(x) = x^2 - a

f'(x) = 2x

これを上の

Xn+1 = Xn - f(Xn)/f'(Xn)

に代入すると

Xn+1 = Xn - (Xn~2 - a)/2Xn

         = Xn/2 + a/2Xn

これを繰り返していくと近似的にルートを求めることができます。

最後に実際にこの漸化式を10回繰り返し近似的にルートを求めるプログラムを紹介したいと思います。今回は10回漸化式を繰り返す方法をとりましたが、自分で近似の精度を決めたい場合はwhile(True)で回し続けて精度がよくなったらbreakで抜けるという方法を取れば良いでしょう。

C言語版

double my_sqrt(double x){
　int i;
　double y,z;
　if(x==0){
　　return 0;
　}
　else{
　　y = 1;
　　for(i=0;i<10;i++){
　　　z = x/y;
　　　y = (y+z)/2;
　　}
　　return y;
　}
}

Java版

public static double my_sqrt(double x){
　if(x==0){
　　return 0;
　}
　else{
　　double y = 1;
　　for(int i=0;i<10;i++){
　　　double z = x/y;
　　　y = (y+z)/2;
　　}
　return y;
　}
}

Python版

def my_sqrt(x):
　if x == 0:
　　return 0
　else:
　　y = 1
　　for i in range(10):
　　　z = x/y
　　　y = (y + z)/2
　　return y

Ruby版

def my_sqrt(x)
　if x == 0
　　return 0
　else
　　y = 1.0
　　for i in 0..9 do
　　　z = x/y
　　　y = (y+z)/2
　　end
　return y
　end
end