(* Staged Matrix Determinant

   Author:  Walid Taha
   Date:    Thu Jun 26 08:23:32 CDT 2003

   Simple implementation of matrix determinant.
   Used to get a symbolic expression of determinant.
*)

let sub (i,j) a (x,y) = 
   a ((if x<i then x else x+1),
      (if y<j then y else y+1));;

let wrap a (x,y) = a.(x).(y);;

let a = Array.make_matrix 3 3 .<0.0>.;;

let lift x = .<x>.;;

for i=0 to 2 do
 for j=0 to 2 do
  a.(i).(j) <- lift (float (if i>j then i+j+5 else (i+1)*(j+j+j+1)))
 done
done;;

let a' = wrap a;;

let rec sum (i,j,f) =
 if i=j then f i
        else .<.~(f i) +. .~(sum (i+1,j,f))>.;;

let det a = 
 let rec det' n a = 
  match n with
    1 -> a (0,0)
  | n -> sum (0,n-1, fun j -> .<.~(lift (if j mod 2 = 0 then 1.0 else -1.0))
                              *. .~(a (0,j))
                              *. .~(det' (n-1) (sub (0,j) a))>.)
 in det' (Array.length a) (wrap a)

let q = Array.make 3 (Array.make 3 .<0.0>.);;

for i=0 to 1 do
 q.(i) <- (Array.make 3 .<0.0>.)
done;;

let det3by3 = .<fun a1 -> fun a2 -> fun a3 -> fun b1 ->
  fun b2 -> fun b3 -> fun c1 -> fun c2 -> fun c3 ->
  .~(q.(0).(0) <- .<a1>.;
     q.(0).(1) <- .<a2>.;
     q.(0).(2) <- .<a3>.;
     q.(1).(0) <- .<b1>.;
     q.(1).(1) <- .<b2>.;
     q.(1).(2) <- .<b3>.;
     q.(2).(0) <- .<c1>.;
     q.(2).(1) <- .<c2>.;
     q.(2).(2) <- .<c3>.;
     det (q))>.;;

(*

det3by3 evaluates to: (* compare to textbook definitions *)

.<fun a1_10 ->
   fun a2_11 ->
    fun a3_12 ->
     fun b1_13 ->
      fun b2_14 ->
       fun b3_15 ->
        fun c1_16 ->
         fun c2_17 ->
          fun c3_18 ->
           (((1 *. a1_10) *.
             (((1 *. b2_14) *. c3_18) +. ((-1 *. b3_15) *. c2_17))) +.
            (((-1 *. a2_11) *.
              (((1 *. b1_13) *. c3_18) +. ((-1 *. b3_15) *. c1_16))) +.
             ((1 *. a3_12) *.
              (((1 *. b1_13) *. c2_17) +. ((-1 *. b2_14) *. c1_16)))))>.
*)