(* rc (read complex) takes a list of zones, each a simplex
   denoted as a list of vertex numbers and builds two lists of
   faces, each sorted by vertex number and segregated by parity.
   This list should be the complex consisting of the boundary. *)


let rc b = List.fold_left
    (fun acc fac -> match fac with (o, vl) ->
   (* Do all faces of this zone with face list fac. *)
   (* Yield a pair of face lists assorted by parity. *)
    (List.fold_left (fun acc1 face ->
      match acc1 with even, odd -> if o
      then (odd, face::even)
      else (face::odd, even))
      acc (sl vl))) ([], []) (List.map pe b);;

let rec cas a b = match a, b with
    | [], [] -> 0 
    | [], a -> -1
    | a, [] ->1
    | p::q, r::s -> let w = compare p r in
       if w = 0 then cas q s else w;;

let ob = List.sort cas;;

let boun q = match (rc q) with a, b -> (ob a, ob b);;


let rec elim q = match q with
   | [], [] -> [], []
   | a, [] -> a, []
   | [], b -> [], b
   | a::b, c::d -> let w = cas a c in if w < 0 then
        match elim (b, c::d) with p, q -> a::p, q else if w> 0 then
        match elim (a::b, d) with p, q -> p, c::q else
        elim (b, d);;

let comb a = elim (boun a);;