ガウスの消去法

　学校の課題でガウスの消去法を実装しました。言語は特に指定されなかったのでOCamlです。

let ( @. ) f g x = f (g x)
let on f g x y = f (g x) (g y)

let ( -^ ) = List.map2 ( -. )
let ( *^ ) x = List.map (( *. ) x)

(* 前進消去 *)
let rec forward_elimination acc matrix =
  (* 枢軸選択 *)
  match List.sort (on compare @@ ( ~-. ) @. abs_float @. List.hd) matrix with
  | [] -> acc
  | (pivot_hd :: pivot_tl) :: rest ->
      let pivot = 1. /. pivot_hd *^ pivot_tl in
      forward_elimination (pivot :: acc)
        @@ List.map (fun (row_hd :: row_tl) -> row_tl -^ row_hd *^ pivot) rest
let forward_elimination = forward_elimination [];;

(* 後退代入 *)
let rec back_substitution acc = function
  | [] -> acc
  | [c] :: rest ->
      back_substitution (c :: acc)
      @@ List.map (fun (c' :: m :: ms) -> (c' -. m *. c) :: ms) rest
let back_substitution = back_substitution [] @. List.map List.rev

let solve = back_substitution @. forward_elimination

(* テスト *)
(* x+y=3, 2x+5=9の解はx=2, y=1 *)
let solve_test1 = solve
    [[1.; 1.; 3.];
     [2.; 5.; 9.]] = [2.; 1.]
let solve_test2 = solve
    [[1.; 2.; 3.];
     [2.; 5.; 5.]] = [5.; -1.]
(* 浮動小数点計算による誤差でテストが通らない *)
let solve_test3 = solve
    [[3.; 1.; 1.];
     [5.; 2.; 1.]] = [1.; -2.]
let solve_test4 = solve
    [[1.; -1.; 1.; -5.];
     [2.; 1.; -3.; 19.];
     [3.; 2.; -1.; 16.]] = [2.; 3.; -4.]

コードを読めば分かる通り、入力が拡大係数行列になっていなければ正しい結果が得られないばかりか実行時エラーが発生してしまいます。
できれば型システムを利用してコンパイル時に弾きたいですね。