let gauss a b =
  let m = Array.length a in
  let n = Array.length a.(0) in

  (* Elementary operation on lines *)
  let exchange i1 i2 =
    if i1 <> i2 then
      let (c,d) = (a.(i1),b.(i1)) in
        a.(i1) <- a.(i2) ;
        a.(i2) <- c ;
        b.(i1) <- b.(i2) ;
        b.(i2) <- d
  in
  let mult i c =
    for j = 0 to n -1 do
      a.(i).(j) <- c *. a.(i).(j)
    done ;
    b.(i) <- c *. b.(i)
  in
  let add i1 i2 c =
    for j = 0 to n - 1 do
      a.(i1).(j) <- a.(i1).(j) +. c *. a.(i2).(j)
    done ;
    b.(i1) <- b.(i1) +. c *. b.(i2)
  in

  let score i =
    (* The score of line i is the first j such that c_ij <> 0. *)
    try
      for j = 0 to n-1 do
        if a.(i).(j) <> 0. then raise (Found j)
      done ; n
    with
      | Found j -> j
  in
  let find i =
    (* We extract the line k>=i with the lowest score. *)
    let best_of (i,si) (k,sk) =
      if si<sk then
        (i,si)
      else
        (k,sk)
    in
    let rec find best k =
      if k = m then best else
        find (best_of best (k,(score k))) (k+1)
    in
      find (i,(score i)) i
  in

    for i = 0 to (min m n) - 1 do
      let (ii,col) = find i in
        if col <> n then
          begin
            (* We'll work with equation ii, pivoting with column col. *)
            exchange i ii ;
            
            (* Normalize *)
            mult i (1./.a.(i).(col)) ;

            (* We want c_{k>i,col} = 0. *)
            for k = i + 1 to m - 1 do
              add k i (-. a.(k).(col)) ;
              a.(k).(col) <- 0.
            done
          end
    done