let gauss a b =
let m = Array.length a in
let n = Array.length a.(0) in
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 =
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 =
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
exchange i ii ;
mult i (1./.a.(i).(col)) ;
for k = i + 1 to m - 1 do
add k i (-. a.(k).(col)) ;
a.(k).(col) <- 0.
done
end
done