大学でニュートン法とDKA法で方程式の解を求める課題が出たので、今回は趣向を変えてStandard MLで実装しました。
実装は以下の通りです。
(* 添字も渡すmap *) fun mapi f = rev o #1 o foldl (fn (x, (result, i)) => (f (i, x) :: result, i + 1)) ([], 0) (* 回数の上限に達するか終了条件を満たすまで計算を反復する *) fun iterateWhile fin f = let fun iterateWhileAux acc 0 x = acc | iterateWhileAux acc n x = if fin x then acc else iterateWhileAux (x :: acc) (n - 1) (f x) in iterateWhileAux [] end local open Complex (* 総乗 *) fun producti f = foldl op* (fromInt 1) o mapi f (* ニュートン法の収束判定 *) fun finNewton eps f x = magnitude (f x) < eps (* DK法の収束判定 *) fun finDK eps f = List.all (finNewton eps f) (* 多項式の微分係数の近似値 *) fun f' c0 j zs zj = producti (fn (i, zi) => if i = j then c0 else zj - zi) zs (* DK法の初期値を求める *) fun aberth (c0 :: c1 :: cs) r = let val n = length (c1 :: cs) (* あまり美しくない *) open Real fun angle k = Math.pi * (2.0 * fromInt k + 0.5) / fromInt n open Complex in List.tabulate (n, fn k => ~c1 / (c0 * fromInt n) + fromReal r * exp (i * fromReal (angle k))) end in (* ホーナー法による多項式の計算 *) fun horner cs z = foldl (fn (c, sum) => sum * z + c) (fromInt 0) cs (* ニュートン法における計算の1ステップ *) fun newton1 f f' x = x - f x / f' x (* 指定された精度で収束するまでニュートン法を反復する *) fun newton eps f f' = iterateWhile (finNewton eps f) (newton1 f f') (* DK法の1ステップ *) fun durandKerner1 (cs as (c0 :: _)) zs = mapi (fn (j, zj) => newton1 (horner cs) (f' c0 j zs) zj) zs (* 指定された精度で収束するまでDK法を反復する *) fun durandKerner eps r cs n = iterateWhile (finDK eps (horner cs)) (durandKerner1 cs) n (aberth cs r) end local open Complex in (* x^2-6+9=0の解をニュートン法で求める *) val q1Newton = newton 1E~6 (horner (map fromInt [1, ~6, 9])) (horner (map fromInt [2, ~6])) 100 (fromReal 0.0) (* x^3-1=0の解をニュートン法で求める *) val q2Newton = map (newton 1E~6 (horner (map fromInt [1, 0, 0, ~1])) (horner (map fromInt [3, 0, 0])) 100) [fromReal 0.5, {real = ~1.0, imag = 1.0}, {real = ~1.0, imag = ~1.0}] (* x^6-4x^4+2x^2-1=0の解をニュートン法で求める *) val q3Newton = map (newton 1E~6 (horner (map fromInt [8, 0, ~4, 0, 2, 0, ~1])) (horner (map fromInt [48, 0, ~16, 0, 4, 0])) 1000) [fromInt 1, fromInt ~1, {real = 1.0, imag = 1.0}, {real = ~1.0, imag = 1.0}, {real = ~1.0, imag = ~1.0}, {real = 1.0, imag = ~1.0}] (* x^2-6+9=0の解をDKA法で求める *) val q1DK = durandKerner 1E~6 1.0 (map fromInt [1, ~6, 9]) 100 (* x^3-1=0の解をDKA法で求める *) val q2DK = durandKerner 1E~6 1.0 (map fromInt [1, 0, 0, ~1]) 100 (* x^6-4x^4+2x^2-1=0の解をDKA法で求める *) val q3DK = durandKerner 1E~6 10.0 (map fromInt [8, 0, ~4, 0, 2, 0, ~1]) 100 end
ここでは複素数計算に適当にググって出てきたhttp://www.cs.unm.edu/~rlpm/code/のライブラリを使っています。 ただし、このライブラリーは除算がバグっていて以下のような修正をしないと正しい値が求められません。
diff -u Complex/Complex.sml Complex_after/Complex.sml --- Complex/Complex.sml 2005-12-17 02:00:43.000000000 +0900 +++ Complex_after/Complex.sml 2014-11-09 21:31:34.000000000 +0900 @@ -88,7 +88,7 @@ fun f x = Real./(x,Real.+(Real.*(a,a),Real.*(b,b))) in {real=f (Real.+(Real.*(x,a),Real.*(y,b))), - imag=f (Real.+(Real.*(y,a),Real.*(x,b)))} + imag=f (Real.-(Real.*(y,a),Real.*(x,b)))} end
このバグのせいで休日が1日潰れてしまいました。 数値計算のプログラムは同じ様な型の値を大量に扱うせいで型の庇護が受けられない上、値を見ただけではどこがおかしいのか分かりづらいのでデバッグは非常につらいものがあります。
Standard MLを選んだばかりに苦労した形になって後味が悪いですが、 OCamlとStandard ML、同じML系の言語にも関わらず構文や命名規則等色々な点が違っていて面白いです。 研究室次第では来年から書く事になるので今から慣れておくのも良さそうですね。