Updated to latest js_of_ocamlHEADmaster

1 files changed, 196 insertions, 222 deletions

diff --git a/script.it/shapes/matrix/Matrix.ml b/script.it/shapes/matrix/Matrix.ml
index 7f1d54b..c71e7c5 100755
--- a/script.it/shapes/matrix/Matrix.ml
+++ b/script.it/shapes/matrix/Matrix.ml
@@ -1,5 +1,4 @@
 open Order
-
 module Order = Order
 
 (*************** Exceptions ***************)
@@ -8,11 +7,8 @@ exception NonSquare
 exception ImproperDimensions
 
 (* Functor so we can Abstract away! *)
-module MakeMatrix (C: EltsI.ORDERED_AND_OPERATIONAL) :
-  (MatrixI.MATRIX with type elt = C.t) =
-struct
-
-
+module MakeMatrix (C : EltsI.ORDERED_AND_OPERATIONAL) :
+  MatrixI.MATRIX with type elt = C.t = struct
   (*************** End Exceptions ***************)
 
   (*************** Types ***************)
@@ -21,7 +17,7 @@ struct
 
   (* A matrix is a pair of dimension (n x p) and a array of arrays
    * the first array is the row (n) and the second the column (p) *)
-  type matrix = (int * int) * (elt array array)
+  type matrix = (int * int) * elt array array
 
   (*************** End Types ***************)
 
@@ -29,12 +25,13 @@ struct
 
   (* catching negative dimensions AND 0 dimensions and too large
    * of a dimension so we don't have to worry about it later *)
-  let empty (rows: int) (columns: int) : matrix =
+  let empty (rows : int) (columns : int) : matrix =
     if rows > 0 && columns > 0 then
       try
-        let m = Array.make_matrix rows columns C.zero in ((rows,columns),m)
-      with _ ->
-        raise ImproperDimensions
+        let m = Array.make_matrix rows columns C.zero in
+        ((rows, columns), m)
+      with
+      | _ -> raise ImproperDimensions
     else (* dimension is negative or 0 *)
       raise ImproperDimensions
 
@@ -50,104 +47,84 @@ struct
    *        |*2 5 6 |
   *)
   (* aside: we don't check whether n < 1 because of our matrix invariant *)
-  let get_row (((n,p),m): matrix) (row: int) : int * elt array =
+  let get_row (((n, p), m) : matrix) (row : int) : int * elt array =
     if row <= n then
       let row' = Array.map (fun x -> x) m.(row - 1) in
       (p, row')
-    else
-      raise (Failure "Row out of bounds.")
+    else raise (Failure "Row out of bounds.")
 
   (* similar to get_row. For m, get_column m 1 will return (2, |1 2|) *)
-  let get_column (((n,p),m): matrix) (column: int) : int * elt array =
-    if column <= p then
-      begin
-        let column' = Array.make n C.zero in
-        for i = 0 to n - 1 do
-          column'.(i) <- m.(i).(column - 1)
-        done;
-        (n, column')
-      end
-    else
-      raise (Failure "Column out of bounds.")
+  let get_column (((n, p), m) : matrix) (column : int) : int * elt array =
+    if column <= p then (
+      let column' = Array.make n C.zero in
+      for i = 0 to n - 1 do
+        column'.(i) <- m.(i).(column - 1)
+      done;
+      (n, column'))
+    else raise (Failure "Column out of bounds.")
 
   (* sets the nth row of the matrix m to the specified array a.
    * This is done IN-PLACE. Therefore the function returns unit.  You should
    * nonetheless enfore immutability whenever possible.  For a clarification on
    * what nth row means, look at comment for get_row above. *)
-  let set_row (((n, p), m): matrix) (row: int) (a: elt array) : unit =
-    if row <= n then
-      begin
-        assert(Array.length a = p);
-        for i = 0 to p - 1 do
-          m.(row - 1).(i) <- a.(i)
-        done;
-      end
-    else
-      raise (Failure "Row out of bounds.")
+  let set_row (((n, p), m) : matrix) (row : int) (a : elt array) : unit =
+    if row <= n then (
+      assert (Array.length a = p);
+      for i = 0 to p - 1 do
+        m.(row - 1).(i) <- a.(i)
+      done)
+    else raise (Failure "Row out of bounds.")
 
   (* Similar to set_row but sets the nth column instead *)
-  let set_column (((n,p),m): matrix) (column: int) (a: elt array) : unit =
-    if column <= p then
-      begin
-        assert(Array.length a = n);
-        for i = 0 to n - 1 do
-          m.(i).(column - 1) <- a.(i)
-        done;
-      end
-    else
-      raise (Failure "Column out of bounds.")
+  let set_column (((n, p), m) : matrix) (column : int) (a : elt array) : unit =
+    if column <= p then (
+      assert (Array.length a = n);
+      for i = 0 to n - 1 do
+        m.(i).(column - 1) <- a.(i)
+      done)
+    else raise (Failure "Column out of bounds.")
 
   (* returns the ij-th element of a matrix (not-zero indexed) *)
-  let get_elt (((n,p),m): matrix) ((i,j): int*int) : elt =
-    if i <= n && j <= p then
-      m.(i - 1).(j - 1)
-    else
-      raise ImproperDimensions
+  let get_elt (((n, p), m) : matrix) ((i, j) : int * int) : elt =
+    if i <= n && j <= p then m.(i - 1).(j - 1) else raise ImproperDimensions
 
   (* Changes the i,j-th element of a matrix to e. Is not zero-indexed, and
    * changes the matrix in place *)
-  let set_elt (((n,p),m): matrix) ((i,j): int*int) (e: elt) : unit =
-    if i <= n && j <= p then
-      m.(i - 1).(j - 1) <- e
-    else
-      raise ImproperDimensions
+  let set_elt (((n, p), m) : matrix) ((i, j) : int * int) (e : elt) : unit =
+    if i <= n && j <= p then m.(i - 1).(j - 1) <- e
+    else raise ImproperDimensions
 
   (* similar to map, but applies to function to the entire matrix
    * Returns a new matrix *)
-  let map (f: elt -> elt) (mat: matrix) : matrix =
-    let (dim,m) = mat in
+  let map (f : elt -> elt) (mat : matrix) : matrix =
+    let dim, m = mat in
     (dim, Array.map (Array.map f) m)
 
   (* Just some wrapping of Array.iter made for Matrices! *)
-  let iter (f: elt -> unit) (mat: matrix) : unit =
+  let iter (f : elt -> unit) (mat : matrix) : unit =
     let _, m = mat in
     Array.iter (Array.iter f) m
 
   (* Just some wrapping of Array.iteri. Useful for pretty
    * printing matrix. The index is (i,j). NOT zero-indexed *)
-  let iteri (f: int -> int -> elt -> unit) (mat: matrix) : unit =
+  let iteri (f : int -> int -> elt -> unit) (mat : matrix) : unit =
     let _, m = mat in
     Array.iteri (fun i row -> Array.iteri (fun j e -> f i j e) row) m
 
   (* folds over each row using base case u and function f *)
   (* could be a bit more efficient? *)
-  let reduce (f: 'a -> elt -> 'a) (u: 'a) (((p,q),m): matrix) : 'a =
+  let reduce (f : 'a -> elt -> 'a) (u : 'a) (((p, q), m) : matrix) : 'a =
     let total = ref u in
     for i = 0 to p - 1 do
       for j = 0 to q - 1 do
-        total := f (!total) m.(i).(j)
-      done;
+        total := f !total m.(i).(j)
+      done
     done;
     !total
 
-  let fold_row ~(f: elt array -> 'b) ((_,m): matrix) : 'b list =
-
-    let call_row acc v = (f v)::acc in
-    Array.fold_left call_row [] m
-    |> List.rev
-
-
-
+  let fold_row ~(f : elt array -> 'b) ((_, m) : matrix) : 'b list =
+    let call_row acc v = f v :: acc in
+    Array.fold_left call_row [] m |> List.rev
 
   (* given two arrays, this will calculate their dot product *)
   (* It seems a little funky, but this is done for efficiency's sake.
@@ -160,87 +137,84 @@ struct
 
      Instead we calculate the length before starting
   *)
-  let dot (v1: elt array) (v2: elt array) : elt =
-    let rec dotting (i: int) (total: elt) : elt =
+  let dot (v1 : elt array) (v2 : elt array) : elt =
+    let rec dotting (i : int) (total : elt) : elt =
       if i = 0 then total
       else
-        let curr = C.multiply v1.(i-1) v2.(i-1) in
-        dotting (i - 1) (C.add curr total) in
-    let len1, len2 = Array.length v1, Array.length v2 in
-    if len1 = len2 then dotting len1 C.zero
-    else raise ImproperDimensions
+        let curr = C.multiply v1.(i - 1) v2.(i - 1) in
+        dotting (i - 1) (C.add curr total)
+    in
+    let len1, len2 = (Array.length v1, Array.length v2) in
+    if len1 = len2 then dotting len1 C.zero else raise ImproperDimensions
 
   (* function to expose the dimensions of a matrix *)
-  let get_dimensions (m: matrix) : (int * int) =
-    let ((x,y), _) = m in (x,y)
+  let get_dimensions (m : matrix) : int * int =
+    let (x, y), _ = m in
+    (x, y)
 
   (*************** End Helper Functions ***************)
 
-
   (*************** Primary Matrix Functions ***************)
 
   (* scales a matrix by the appropriate factor *)
-  let scale (m: matrix) (sc: elt) : matrix = map (C.multiply sc) m
+  let scale (m : matrix) (sc : elt) : matrix = map (C.multiply sc) m
 
   (* Generates a matrix from a list of lists. The inners lists are the rows *)
   let from_list (lsts : elt list list) : matrix =
-    let check_length (length: int) (lst: elt list) : int =
-      if List.length lst = length then length
-      else raise ImproperDimensions in
+    let check_length (length : int) (lst : elt list) : int =
+      if List.length lst = length then length else raise ImproperDimensions
+    in
     let p = List.length lsts in
     match lsts with
     | [] -> raise ImproperDimensions
-    | hd::tl ->
-      let len = List.length hd in
-      if List.fold_left check_length len tl = len then
-        ((p,len),Array.map Array.of_list (Array.of_list lsts))
-      else
-        raise ImproperDimensions
+    | hd :: tl ->
+        let len = List.length hd in
+        if List.fold_left check_length len tl = len then
+          ((p, len), Array.map Array.of_list (Array.of_list lsts))
+        else raise ImproperDimensions
 
   (* Generates a matrix from a list of lists. The inners lists are the rows *)
   let from_array (arrs : elt array array) : matrix =
-    let check_length (length: int) (arr: elt array) : unit =
-      if Array.length arr = length then ()
-      else raise ImproperDimensions in
+    let check_length (length : int) (arr : elt array) : unit =
+      if Array.length arr = length then () else raise ImproperDimensions
+    in
     let p = Array.length arrs in
     match Array.length arrs with
     | 0 -> raise ImproperDimensions
     | _ ->
-      let len = Array.length (Array.get arrs 0) in
-      Array.iter (check_length len) arrs;
-      ((p, len), arrs)
+        let len = Array.length (Array.get arrs 0) in
+        Array.iter (check_length len) arrs;
+        ((p, len), arrs)
 
   (* Adds two matrices. They must have the same dimensions *)
-  let add ((dim1,m1): matrix) ((dim2,m2): matrix) : matrix =
-    if dim1 = dim2 then
+  let add ((dim1, m1) : matrix) ((dim2, m2) : matrix) : matrix =
+    if dim1 = dim2 then (
       let n, p = dim1 in
-      let (dim', sum_m) = empty n p in
+      let dim', sum_m = empty n p in
       for i = 0 to n - 1 do
         for j = 0 to p - 1 do
           sum_m.(i).(j) <- C.add m1.(i).(j) m2.(i).(j)
-        done;
+        done
       done;
-      (dim',sum_m)
-    else
-      raise ImproperDimensions
-
+      (dim', sum_m))
+    else raise ImproperDimensions
 
   (* Multiplies two matrices. If the matrices have dimensions m x n and p x q, n
    * and p must be equal, and the resulting matrix will have dimension n x q *)
-  let mult (matrix1: matrix) (matrix2: matrix) : matrix =
-    let ((m,n), _), ((p,q), _) = matrix1, matrix2 in
-    if n = p then
-      let (dim, result) = empty m q in
+  let mult (matrix1 : matrix) (matrix2 : matrix) : matrix =
+    let ((m, n), _), ((p, q), _) = (matrix1, matrix2) in
+    if n = p then (
+      let dim, result = empty m q in
       for i = 0 to m - 1 do
         for j = 0 to q - 1 do
-          let (_,row), (_,column) = get_row matrix1 (i + 1),
-                                    get_column matrix2 (j + 1) in
+          let (_, row), (_, column) =
+            (get_row matrix1 (i + 1), get_column matrix2 (j + 1))
+          in
           result.(i).(j) <- dot row column
-        done;
+        done
       done;
-      (dim,result)
-    else
-      raise ImproperDimensions
+      (dim, result))
+    else raise ImproperDimensions
 
   (*************** Helper Functions for Row Reduce ***************)
 
@@ -294,7 +268,8 @@ struct
 
   (* Compares two elements in an elt array and returns the greater and its
    * index. Is a helper function for find_max_col_index *)
-  let compare_helper (e1: elt) (e2: elt) (ind1: int) (ind2: int) : (elt*int) =
+  let compare_helper (e1 : elt) (e2 : elt) (ind1 : int) (ind2 : int) : elt * int
+      =
     match C.compare e1 e2 with
     | Equal -> (e2, ind2)
     | Greater -> (e1, ind1)
@@ -303,50 +278,53 @@ struct
   (* Finds the element with the greatest absolute value in a column. Is not
    * 0-indexed. If two elements are both the maximum value, returns the one with
    * the lowest index. Returns None if this element is zero (if column is all 0)
-  *)
-  let find_max_col_index (array1: elt array) (start_index: int) : int option =
-    let rec find_index (max_index: int) (curr_max: elt) (curr_index: int)
-        (arr: elt array) =
+   *)
+  let find_max_col_index (array1 : elt array) (start_index : int) : int option =
+    let rec find_index (max_index : int) (curr_max : elt) (curr_index : int)
+        (arr : elt array) =
       if curr_index = Array.length arr then
-        (if curr_max = C.zero then None
-         else Some (max_index+1)) (* Arrays are 0-indexed but matrices aren't *)
+        if curr_max = C.zero then None
+        else Some (max_index + 1) (* Arrays are 0-indexed but matrices aren't *)
       else
-        (match C.compare arr.(curr_index) C.zero with
-         | Equal -> find_index max_index curr_max (curr_index+1) arr
-         | Greater ->
-           (let (el, index) = compare_helper (arr.(curr_index))
-                curr_max curr_index max_index in
-            find_index index el (curr_index+1) arr)
-         | Less ->
-           (let abs_curr_elt = C.subtract C.zero arr.(curr_index) in
-            let (el, index) = compare_helper abs_curr_elt curr_max curr_index
-                max_index in
-            find_index index el (curr_index+1) arr))
+        match C.compare arr.(curr_index) C.zero with
+        | Equal -> find_index max_index curr_max (curr_index + 1) arr
+        | Greater ->
+            let el, index =
+              compare_helper arr.(curr_index) curr_max curr_index max_index
+            in
+            find_index index el (curr_index + 1) arr
+        | Less ->
+            let abs_curr_elt = C.subtract C.zero arr.(curr_index) in
+            let el, index =
+              compare_helper abs_curr_elt curr_max curr_index max_index
+            in
+            find_index index el (curr_index + 1) arr
     in
-    find_index 0 C.zero (start_index -1) array1
+    find_index 0 C.zero (start_index - 1) array1
 
   (* Basic row operations *)
   (* Scales a row by sc *)
-  let scale_row (m: matrix) (num: int) (sc: elt) : unit =
-    let (_, row) = get_row m num in
+  let scale_row (m : matrix) (num : int) (sc : elt) : unit =
+    let _, row = get_row m num in
     let new_row = Array.map (C.multiply sc) row in
     set_row m num new_row
 
   (* Swaps two rows of a matrix *)
-  let swap_row (m: matrix) (r1: int) (r2: int) : unit =
-    let (len1, row1) = get_row m r1 in
-    let (len2, row2) = get_row m r2 in
+  let swap_row (m : matrix) (r1 : int) (r2 : int) : unit =
+    let len1, row1 = get_row m r1 in
+    let len2, row2 = get_row m r2 in
     let _ = assert (len1 = len2) in
     let _ = set_row m r1 row2 in
     let _ = set_row m r2 row1 in
     ()
 
   (* Subtracts a multiple of r2 from r1 *)
-  let sub_mult (m: matrix) (r1: int) (r2: int) (sc: elt) : unit =
-    let (len1, row1) = get_row m r1 in
-    let (len2, row2) = get_row m r2 in
+  let sub_mult (m : matrix) (r1 : int) (r2 : int) (sc : elt) : unit =
+    let len1, row1 = get_row m r1 in
+    let len2, row2 = get_row m r2 in
     let _ = assert (len1 = len2) in
-    for i = 0 to len1 - 1 do (* Arrays are 0-indexed *)
+    for i = 0 to len1 - 1 do
+      (* Arrays are 0-indexed *)
       row1.(i) <- C.subtract row1.(i) (C.multiply sc row2.(i))
     done;
     set_row m r1 row1
@@ -355,175 +333,171 @@ struct
 
   (* Returns the row reduced form of a matrix. Is not done in place, but creates
    * a new matrix *)
-  let row_reduce (mat: matrix) : matrix =
-    let[@tailcall] rec row_reduce_h (n_row: int) (n_col: int) (mat2: matrix) : unit =
-      let ((num_row, _), _) = mat2 in
-      if (n_col = num_row + 1) then ()
+  let row_reduce (mat : matrix) : matrix =
+    let rec row_reduce_h (n_row : int) (n_col : int) (mat2 : matrix) : unit =
+      let (num_row, _), _ = mat2 in
+      if n_col = num_row + 1 then ()
       else
-        let (_,col) = get_column mat2 n_col in
+        let _, col = get_column mat2 n_col in
         match find_max_col_index col n_row with
-        | None (* Column all 0s *) -> row_reduce_h n_row (n_col+1) mat2
+        | None (* Column all 0s *) ->
+            (row_reduce_h [@tailcall]) n_row (n_col + 1) mat2
         | Some index ->
-          begin
             swap_row mat2 index n_row;
             let pivot = get_elt mat2 (n_row, n_col) in
-            scale_row mat2 (n_row) (C.divide C.one pivot);
+            scale_row mat2 n_row (C.divide C.one pivot);
             for i = 1 to num_row do
-              if i <> n_row then sub_mult mat2 i n_row (get_elt mat2 (i,n_col))
+              if i <> n_row then sub_mult mat2 i n_row (get_elt mat2 (i, n_col))
             done;
-            row_reduce_h (n_row+1) (n_col+1) mat2
-          end
+            (row_reduce_h [@tailcall]) (n_row + 1) (n_col + 1) mat2
     in
     (* Copies the matrix *)
-    let ((n,p),m) = mat in
-    let (dim,mat_cp) = empty n p in
+    let (n, p), m = mat in
+    let dim, mat_cp = empty n p in
     for i = 0 to n - 1 do
       for j = 0 to p - 1 do
         mat_cp.(i).(j) <- m.(i).(j)
-      done;
+      done
     done;
-    let _ = row_reduce_h 1 1 (dim,mat_cp) in (dim,mat_cp)
+    let _ = row_reduce_h 1 1 (dim, mat_cp) in
+    (dim, mat_cp)
 
   (*************** End Main Functions ***************)
 
   (*************** Optional module functions ***************)
 
   (* calculates the trace of a matrix *)
-  let trace (((n,p),m): matrix) : elt =
-    let rec build (elt: elt) (i: int) =
-      if i > -1 then
-        build (C.add m.(i).(i) elt) (i - 1)
-      else
-        elt in
-    if n = p then build C.zero (n - 1)
-    else raise ImproperDimensions
+  let trace (((n, p), m) : matrix) : elt =
+    let rec build (elt : elt) (i : int) =
+      if i > -1 then build (C.add m.(i).(i) elt) (i - 1) else elt
+    in
+    if n = p then build C.zero (n - 1) else raise ImproperDimensions
 
   (* calculates the transpose of a matrix and retuns a new one *)
-  let transpose (((n,p),m): matrix) =
-    let (dim,m') = empty p n in
+  let transpose (((n, p), m) : matrix) =
+    let dim, m' = empty p n in
     for i = 0 to n - 1 do
       for j = 0 to p - 1 do
         m'.(j).(i) <- m.(i).(j)
-      done;
+      done
     done;
-    assert(dim = (p,n));
-    ((p,n),m')
+    assert (dim = (p, n));
+    ((p, n), m')
 
   (* Returns the inverse of a matrix. Uses a pretty simple algorithm *)
-  let inverse (mat: matrix) : matrix =
-    let ((n, p), _) = mat in
-    if n = p then
+  let inverse (mat : matrix) : matrix =
+    let (n, p), _ = mat in
+    if n = p then (
       (* create augmented matrix *)
-      let augmented = empty n (2*n) in
+      let augmented = empty n (2 * n) in
       for i = 1 to n do
-        let (dim,col) = get_column mat i in
+        let dim, col = get_column mat i in
         let arr = Array.make n C.zero in
-        begin
-          assert(dim = n);
-          arr.(i-1) <- C.one;
-          set_column augmented i col;
-          set_column augmented (n + i) arr
-        end
+        assert (dim = n);
+        arr.(i - 1) <- C.one;
+        set_column augmented i col;
+        set_column augmented (n + i) arr
       done;
       let augmented' = row_reduce augmented in
       (* create the inverted matrix and fill in with appropriate values *)
       let inverse = empty n n in
       for i = 1 to n do
-        let (dim, col) = get_column augmented' (n + i) in
-        let _ = assert(dim = n) in
+        let dim, col = get_column augmented' (n + i) in
+        let _ = assert (dim = n) in
         let _ = set_column inverse i col in
         ()
       done;
-      inverse
-    else
-      raise NonSquare
+      inverse)
+    else raise NonSquare
 
   (***************** HELPER FUNCTIONS FOR DETERMINANT *****************)
   (* creates an identity matrix of size n*)
-  let create_identity (n:int) : matrix =
-    let (dim,m) = empty n n in
+  let create_identity (n : int) : matrix =
+    let dim, m = empty n n in
     for i = 0 to n - 1 do
       m.(i).(i) <- C.one
     done;
-    (dim,m)
+    (dim, m)
 
   (* Finds the index of the maximum value of an array *)
-  let find_max_index (arr: elt array) (start_index : int) : int =
-    let rec find_index (max_index: int) (curr_index: int) =
-      if curr_index = Array.length arr then max_index+1
+  let find_max_index (arr : elt array) (start_index : int) : int =
+    let rec find_index (max_index : int) (curr_index : int) =
+      if curr_index = Array.length arr then max_index + 1
       else
         match C.compare arr.(curr_index) arr.(max_index) with
         | Equal | Less -> find_index max_index (curr_index + 1)
-        | Greater -> find_index curr_index (curr_index + 1) in
+        | Greater -> find_index curr_index (curr_index + 1)
+    in
     find_index (start_index - 1) start_index
 
   (* Creates the pivoting matrix for A. Returns swqps. Adapted from
    * http://rosettacode.org/wiki/LU_decomposition#Common_Lisp *)
-  let pivotize (((n,p),m): matrix) : matrix * int =
-    if n = p then
+  let pivotize (((n, p), m) : matrix) : matrix * int =
+    if n = p then (
       let swaps = ref 0 in
       let pivot_mat = create_identity n in
       for j = 1 to n do
-        let (_,col) = get_column ((n,p),m) j in
+        let _, col = get_column ((n, p), m) j in
         let max_index = find_max_index col j in
-        if max_index <> j then
-          (swaps := !swaps + 1; swap_row pivot_mat max_index j)
+        if max_index <> j then (
+          swaps := !swaps + 1;
+          swap_row pivot_mat max_index j)
         else ()
       done;
-      (pivot_mat,!swaps)
+      (pivot_mat, !swaps))
     else raise ImproperDimensions
 
   (* decomposes a matrix into a lower triangualar, upper triangualar
    * and a pivot matrix. It returns (L,U,P). Adapted from
    * http://rosettacode.org/wiki/LU_decomposition#Common_Lisp *)
-  let lu_decomposition (((n,p),m): matrix) : (matrix*matrix*matrix)*int =
-    if n = p then
-      let mat = ((n,p),m) in
-      let lower, upper, (pivot,s) = empty n n, empty n n, pivotize mat in
-      let (_ ,l),(_ ,u), _ = lower,upper,pivot in
-      let ((_, _),mat') = mult pivot mat in
+  let lu_decomposition (((n, p), m) : matrix) : (matrix * matrix * matrix) * int
+      =
+    if n = p then (
+      let mat = ((n, p), m) in
+      let lower, upper, (pivot, s) = (empty n n, empty n n, pivotize mat) in
+      let (_, l), (_, u), _ = (lower, upper, pivot) in
+      let (_, _), mat' = mult pivot mat in
       for j = 0 to n - 1 do
         l.(j).(j) <- C.one;
         for i = 0 to j do
           let sum = ref C.zero in
           for k = 0 to i - 1 do
-            sum := C.add (!sum) (C.multiply u.(k).(j) l.(i).(k))
+            sum := C.add !sum (C.multiply u.(k).(j) l.(i).(k))
           done;
-          u.(i).(j) <- C.subtract mat'.(i).(j) (!sum)
+          u.(i).(j) <- C.subtract mat'.(i).(j) !sum
         done;
         for i = j to n - 1 do
           let sum = ref C.zero in
           for k = 0 to j - 1 do
-            sum := C.add (!sum) (C.multiply u.(k).(j) l.(i).(k))
+            sum := C.add !sum (C.multiply u.(k).(j) l.(i).(k))
           done;
-          let sub = C.subtract mat'.(i).(j) (!sum) in
+          let sub = C.subtract mat'.(i).(j) !sum in
           l.(i).(j) <- C.divide sub u.(j).(j)
-        done;
+        done
       done;
-      (lower,upper,pivot),s
+      ((lower, upper, pivot), s))
     else raise ImproperDimensions
 
   (* Computes the determinant of a matrix *)
-  let determinant (m: matrix) : elt =
+  let determinant (m : matrix) : elt =
     try
-      let ((n,p), _) = m in
+      let (n, p), _ = m in
       if n = p then
-        let rec triangualar_det (a,mat) curr_index acc =
+        let rec triangualar_det (a, mat) curr_index acc =
           if curr_index < n then
             let acc' = C.multiply mat.(curr_index).(curr_index) acc in
-            triangualar_det (a,mat) (curr_index + 1) acc'
-          else acc in
-        let ((dim1,l),(dim2,u), _),s = lu_decomposition m in
-        let det1, det2 = triangualar_det (dim1,l) 0 C.one,
-                         triangualar_det (dim2,u) 0 C.one in
+            triangualar_det (a, mat) (curr_index + 1) acc'
+          else acc
+        in
+        let ((dim1, l), (dim2, u), _), s = lu_decomposition m in
+        let det1, det2 =
+          (triangualar_det (dim1, l) 0 C.one, triangualar_det (dim2, u) 0 C.one)
+        in
         if s mod 2 = 0 then C.multiply det1 det2
         else C.subtract C.zero (C.multiply det1 det2)
       else raise ImproperDimensions
     with
     | _ -> C.zero
 
-
   (*************** Optional module functions ***************)
-
-
 end