任意精度の有理数計算によるNewton-Raphson法
# 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"