読者です 読者をやめる 読者になる 読者になる

任意精度の有理数計算によるNewton-Raphson法

OCaml
# open Num;;
# (* f'(x)を求める *)
  let deriv f =
    let dx = Int 10 **/ Int (-50) in
      fun x -> (f (x +/ dx) -/ f x) // dx;;
val deriv : (Num.num -> Num.num) -> Num.num -> Num.num = <fun>

# (* f(x), f(f(x)), f(f(f(x))), ..., を計算していき、ステップ間の差が閾値より小さくなったときの値を返す *)
  let fixpoint f init =
    let threshold = Int 10 **/ Int (-50) in
    let rec loop x =
      let next = f x in
      if abs_num (x -/ next) </ threshold then x
      else loop next
    in loop init;;
val fixpoint : (Num.num -> Num.num) -> Num.num -> Num.num = <fun>
# (* ある点における接線と、x軸との交点のx座標 (x - f(x) / f'(x)) を返す *) 
  let newton_transform f = fun x -> x -/ f x // (deriv f x);;
val newton_transform : (Num.num -> Num.num) -> Num.num -> Num.num = <fun>
# (* 接線とx軸との交点が動かなくなったらおわり *)
  let newton_method f guess = fixpoint (newton_transform f) guess;;
val newton_method : (Num.num -> Num.num) -> Num.num -> Num.num = <fun>

平方根を求めてみる。

# let square_root a = newton_method (fun x -> x */ x -/ a) (Int 1);;
val square_root : Num.num -> Num.num = <fun>
# let s5 = square_root (Int 5);;
val s5 : Num.num = Ratio <abstr>
# approx_num_fix 50 s5;;
- : string = "+2.23606797749978969640917366873127623544061835961153"
# approx_num_fix 50 (s5 */ s5);;
- : string = "+5.00000000000000000000000000000000000000000000000000"