(* ds2's/2 of a random simplex. *)
let pt px = let n = (Array.length px) - 1 in
let g ds2 = Array.init n (fun i -> Array.init n (fun j ->
   ds2.(n).(i) +. ds2.(n).(j) -. ds2.(i).(j))) in 
let ds2a = Array.init (n+1) (fun i -> Array.init (n+1) (fun j ->
    if i < j then (Random.float 1.) else 0.)) in
for i = 0 to n do for j = 0 to i-1 do
    ds2a.(i).(j) <- ds2a.(j).(i) done done;
let pa a =
let pma p a = Array.init (Array.length p) (fun j -> a.(p.(j))) in
pma px (Array.map (pma px) a) in
let ds2b = pa ds2a in
let ga = g ds2a and gb = g ds2b in
let gia = Matrix.inv ga and gib = Matrix.inv gb
and aug a = let n = Array.length a in let lr = Array.make (n+1) 0. in
        Array.init (n+1) (fun i -> if i < n then
        Array.init (n+1) (fun j -> let e = (if j < n then a.(i).(j) else
            Array.fold_left (-.) 0. a.(i)) in lr.(j) <- lr.(j) -. e; e) else lr) in
let agia = aug gia and agib = aug gib in
let pagia = pa agia in
let pm nm a = let pa a = for i = 0 to (Array.length a)-1 do
    print_float a.(i); print_char ' ' done;
    print_newline () in
  print_string nm; print_newline ();
    for i = 0 to (Array.length a)-1 do pa a.(i) done in

pm "ds2a" ds2a; pm "ds2b" ds2b;
pm "ga" ga; pm "gb" gb;
pm "gia" gia; pm "gib" gib;
pm "agia" agia; pm "pagia" pagia;
pm "agib" agib;
print_float (let m = ref 0. in for i = 0 to n do for j = 0 to n do
    m := max !m (abs_float (agib.(i).(j) -. pagia.(i).(j))) done done;
    !m);;

aug [|[|3.; 5.|]; [|5.; 4.|]|];;

pt [|2; 3; 4; 1; 0|];;
pt [|2; 1; 0; 4; 3|];;
pt [|2; 3; 0; 1; 4|];;
pt [|0; 1; 2; 3; 4|];;