diff --git a/.merlin b/.merlin
deleted file mode 100644
index 0f3f94b..0000000
--- a/.merlin
+++ /dev/null
@@ -1,9 +0,0 @@
-S src
-S lib
-S dgraph
-B src
-B lib
-B dgraph
-B .
-PKG lablgtk2
-FLG -w +a -w -4 -w -44 -w -48 -w -29
diff --git a/CHANGES.md b/CHANGES.md
index cd8beac..8e837f3 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,8 +1,18 @@
 
+  - [Classic]: new functions [cycle] and [grid]
+  - [Eulerian]: Eulerian paths (new module)
+  - [Components]: strong articulation points (see functors [Connectivity]
+    and [BiConnectivity]) (Timothy Bourke)
+  - [Dominator]: non-trivial dominators (Timothy Bourke)
+  - #31: fixed documentation of [map_vertex]: the supplied function
+    must be injective
+  - #110: ensure that map_vertex applies the function only once per vertex
+
 # 2.0.0 (October 2, 2020)
 
   - port to dune and opam 2.0
-  - :exclamation: opam package now split into two packages: ocamlgraph and ocamlgraph_gtk
+  - :exclamation: opam package now split into two packages: ocamlgraph
+    and ocamlgraph_gtk
   - [WeakTopological] fixed incorrect use of generic hash tables
     (#99, Tomáš Dacík)
   - [Oper] fixed transitive_reduction (#91)
@@ -12,7 +22,11 @@
   - :exclamation: [Coloring] functions now fail if the graph is directed
   - :exclamation: [Coloring] now uses a single, global exception [NoColoring]
   - [Coloring] new function two_color to 2-color a graph (or fail)
-  - :exclamation: [Fixpoint] Take initial labeling of nodes into account (Johannes Kloos)
+  - :exclamation: [Fixpoint] Take initial labeling of nodes into
+    account (Johannes Kloos)
+  - :exclamation: [Dominator.Make_graph] now accepts a signature that
+    is Builder-compatible
+
 
 # 1.8.8, October 17, 2017
 
diff --git a/CREDITS b/CREDITS
index 24ee8c4..a7fa5f4 100644
--- a/CREDITS
+++ b/CREDITS
@@ -2,8 +2,8 @@
 Contributors, in alphabetic order 
 
   Sylvain Conchon
-  Jean-Christophe Filli�tre
-  Fran�ois Pottier
+  Jean-Christophe Filliâtre
+  François Pottier
   Julien Signoles
   Vincent Simonet
   Matthieu Sozeau
diff --git a/debian/changelog b/debian/changelog
index 3fd268a..9412550 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,10 @@
+ocamlgraph (186+git20220421.1.ef366e7-1) UNRELEASED; urgency=low
+
+  * New upstream snapshot.
+  * Drop patch encoding-utf8, present upstream.
+
+ -- Debian Janitor <janitor@jelmer.uk>  Thu, 08 Sep 2022 17:18:00 -0000
+
 ocamlgraph (2.0.0-3) unstable; urgency=medium
 
   * Team upload
diff --git a/debian/patches/0002-avoid-referencing-non-existent-icon-ed_icon.xpm.patch b/debian/patches/0002-avoid-referencing-non-existent-icon-ed_icon.xpm.patch
index 43de7c0..ae3921b 100644
--- a/debian/patches/0002-avoid-referencing-non-existent-icon-ed_icon.xpm.patch
+++ b/debian/patches/0002-avoid-referencing-non-existent-icon-ed_icon.xpm.patch
@@ -6,10 +6,10 @@ Subject: avoid referencing non existent icon ed_icon.xpm
  editor/ed_main.ml | 4 +++-
  1 file changed, 3 insertions(+), 1 deletion(-)
 
-diff --git a/editor/ed_main.ml b/editor/ed_main.ml
-index b7309cd..b48ef1a 100644
---- a/editor/ed_main.ml
-+++ b/editor/ed_main.ml
+Index: ocamlgraph/editor/ed_main.ml
+===================================================================
+--- ocamlgraph.orig/editor/ed_main.ml
++++ ocamlgraph/editor/ed_main.ml
 @@ -895,7 +895,7 @@ described in file LICENSE.
  This software is distributed in the hope that it will be useful, 
  but WITHOUT ANY WARRANTY; without even the implied warranty of 
diff --git a/debian/patches/encoding-utf8 b/debian/patches/encoding-utf8
deleted file mode 100644
index 9fc67e9..0000000
--- a/debian/patches/encoding-utf8
+++ /dev/null
@@ -1,15 +0,0 @@
-Index: ocamlgraph/CREDITS
-===================================================================
---- ocamlgraph.orig/CREDITS	2020-12-29 18:20:47.866850565 +0100
-+++ ocamlgraph/CREDITS	2020-12-29 18:20:47.000000000 +0100
-@@ -2,8 +2,8 @@
- Contributors, in alphabetic order 
- 
-   Sylvain Conchon
--  Jean-Christophe Filli�tre
--  Fran�ois Pottier
-+  Jean-Christophe Filliâtre
-+  François Pottier
-   Julien Signoles
-   Vincent Simonet
-   Matthieu Sozeau
diff --git a/debian/patches/series b/debian/patches/series
index 2899244..4a6e056 100644
--- a/debian/patches/series
+++ b/debian/patches/series
@@ -1,2 +1 @@
-encoding-utf8
 0002-avoid-referencing-non-existent-icon-ed_icon.xpm.patch
diff --git a/ocamlgraph.opam b/ocamlgraph.opam
index 375499d..360f530 100644
--- a/ocamlgraph.opam
+++ b/ocamlgraph.opam
@@ -15,11 +15,13 @@ tags: [
   "imperative"
 ]
 homepage: "https://github.com/backtracking/ocamlgraph/"
+doc: "https://backtracking.github.io/ocamlgraph"
 bug-reports: "https://github.com/backtracking/ocamlgraph/issues/new"
 depends: [
-  "ocaml"
+  "ocaml" {>= "4.03.0"}
   "stdlib-shims"
   "dune" {>= "2.0"}
+  "graphics" {with-test}
 ]
 build: [
   ["dune" "subst"] {pinned}
diff --git a/ocamlgraph_gtk.opam b/ocamlgraph_gtk.opam
index 6a56a18..07998e5 100644
--- a/ocamlgraph_gtk.opam
+++ b/ocamlgraph_gtk.opam
@@ -15,14 +15,16 @@ tags: [
   "imperative"
 ]
 homepage: "https://github.com/backtracking/ocamlgraph/"
+doc: "https://backtracking.github.io/ocamlgraph"
 bug-reports: "https://github.com/backtracking/ocamlgraph/issues/new"
 depends: [
-  "ocaml"
+  "ocaml" {>= "4.03.0"}
   "stdlib-shims"
   "lablgtk"
   "conf-gnomecanvas"
-  "ocamlgraph"
+  "ocamlgraph" {>= "2.0.0"}
   "dune" {>= "2.0"}
+  "graphics" {with-test}
 ]
 build: [
   ["dune" "subst"] {pinned}
diff --git a/src/blocks.ml b/src/blocks.ml
index 1f9331f..1436c87 100644
--- a/src/blocks.ml
+++ b/src/blocks.ml
@@ -235,8 +235,10 @@ module Unlabeled(V: COMPARABLE)(HM: HM with type key = V.t) = struct
   let succ g v = S.elements (HM.find_and_raise v g "[ocamlgraph] succ")
   let succ_e g v = fold_succ_e (fun e l -> e :: l) g v []
 
-  let map_vertex f =
-    HM.map (fun v s -> f v, S.fold (fun v s -> S.add (f v) s) s S.empty)
+  let map_vertex f g =
+    let module MV = Util.Memo(V) in
+    let f = MV.memo f in
+    HM.map (fun v s -> f v, S.fold (fun v s -> S.add (f v) s) s S.empty) g
 
   module I = struct
     type t = S.t HM.t
@@ -348,9 +350,11 @@ struct
   let succ g v = fold_succ (fun w l -> w :: l) g v []
   let succ_e g v = fold_succ_e (fun e l -> e :: l) g v []
 
-  let map_vertex f =
+  let map_vertex f g =
+    let module MV = Util.Memo(V) in
+    let f = MV.memo f in
     HM.map
-      (fun v s -> f v, S.fold (fun (v, l) s -> S.add (f v, l) s) s S.empty)
+      (fun v s -> f v, S.fold (fun (v, l) s -> S.add (f v, l) s) s S.empty) g
 
   module I = struct
     type t = S.t HM.t
@@ -561,12 +565,15 @@ module BidirectionalUnlabeled(V:COMPARABLE)(HM:HM with type key = V.t) = struct
   let succ g v = S.elements (snd (HM.find_and_raise v g "[ocamlgraph] succ"))
   let succ_e g v = fold_succ_e (fun e l -> e :: l) g v []
 
-  let map_vertex f =
+  let map_vertex f g =
+    let module MV = Util.Memo(V) in
+    let f = MV.memo f in
     HM.map
       (fun v (s1,s2) ->
          f v,
          (S.fold (fun v s -> S.add (f v) s) s1 S.empty,
           S.fold (fun v s -> S.add (f v) s) s2 S.empty))
+      g
 
   module I = struct
     (* we keep sets for both incoming and outgoing edges *)
@@ -703,12 +710,15 @@ struct
   let succ g v = fold_succ (fun w l -> w :: l) g v []
   let succ_e g v = fold_succ_e (fun e l -> e :: l) g v []
 
-  let map_vertex f =
+  let map_vertex f g =
+    let module MV = Util.Memo(V) in
+    let f = MV.memo f in
     HM.map
       (fun v (s1,s2) ->
          f v,
          (S.fold (fun (v, l) s -> S.add (f v, l) s) s1 S.empty,
           S.fold (fun (v, l) s -> S.add (f v, l) s) s2 S.empty))
+    g
 
   module I = struct
     type t = (S.t * S.t) HM.t
diff --git a/src/classic.ml b/src/classic.ml
index 2d66756..9350459 100644
--- a/src/classic.ml
+++ b/src/classic.ml
@@ -19,16 +19,21 @@
 
 module type S = sig
   type graph
+  type vertex
   val divisors : int -> graph
   val de_bruijn : int -> graph
   val vertex_only : int -> graph
   val full : ?self:bool -> int -> graph
+  val cycle : int -> graph * vertex array
+  val grid : n:int -> m:int -> graph * vertex array array
 end
 
 module Generic(B : Builder.INT) = struct
 
   type graph = B.G.t
 
+  type vertex = B.G.V.t
+
   let divisors n =
     if n < 2 then invalid_arg "divisors";
     let v = Array.init (n + 1) (fun i -> B.G.V.create i) in
@@ -78,6 +83,28 @@ module Generic(B : Builder.INT) = struct
            g)
       (fold_for 1 n (fun g i -> B.add_vertex g v.(i)) (B.empty ()))
 
+  let cycle n =
+    if n < 0 then invalid_arg "cycle";
+    let v = Array.init n (fun i -> B.G.V.create i) in
+    let g = Array.fold_left B.add_vertex (B.empty ()) v in
+    let rec loop g i =
+      if i = n then g
+      else let g = B.add_edge g v.(i) v.((i+1) mod n) in loop g (i+1) in
+    loop g 0, v
+
+  let grid ~n ~m =
+    if n < 0 || m < 0 then invalid_arg "grid";
+    let create i j = B.G.V.create (m * i + j) in
+    let v = Array.init n (fun i -> Array.init m (fun j -> create i j)) in
+    let g = Array.fold_left (Array.fold_left B.add_vertex) (B.empty ()) v in
+    let rec loop g i j =
+      if i = n then g
+      else if j = m then loop g (i+1) 0
+      else let g = if j < m-1 then B.add_edge g v.(i).(j) v.(i).(j+1) else g in
+           let g = if i < n-1 then B.add_edge g v.(i).(j) v.(i+1).(j) else g in
+           loop g i (j+1) in
+    loop g 0 0, v
+
 end
 
 module P (G : Sig.P with type V.label = int) = Generic(Builder.P(G))
diff --git a/src/classic.mli b/src/classic.mli
index 0e85ab0..a21e9c0 100644
--- a/src/classic.mli
+++ b/src/classic.mli
@@ -15,14 +15,14 @@
 (*                                                                        *)
 (**************************************************************************)
 
-(* $Id: classic.mli,v 1.12 2005-02-25 13:54:33 signoles Exp $ *)
-
 (** Some classic graphs *)
 
 module type S = sig
 
   type graph
 
+  type vertex
+
   val divisors : int -> graph
   (** [divisors n] builds the graph of divisors.
       Vertices are integers from [2] to [n]. [i] is connected to [j] if
@@ -45,10 +45,24 @@ module type S = sig
       The optional argument [self] indicates if loop edges should be added
       (default value is [true]). *)
 
+  val cycle : int -> graph * vertex array
+  (** [cycle n] builds a graph that is a cycle with [n] vertices.
+      Vertices are labelled with [0,1,...,n-1] and there is an edge from
+      vertex [i] to vertex [(i+1) mod n].
+      Vertices are also returned in an array for convenience. *)
+
+  val grid : n:int -> m:int -> graph * vertex array array
+  (** [grid n m] builds a grid graph with [n*m] vertices, with edges
+      from vertex [(i,j)] to vertices [(i+1,j)] and [(i,j+1)] (and no
+      wrapping around). Vertex [(i,j)] is labelled with [i*m+j].
+      Vertices are also returned in a [n*m] matrix for convenience. *)
+
 end
 
-module P (G : Sig.P with type V.label = int) : S with type graph = G.t
+module P (G : Sig.P with type V.label = int) :
+   S with type graph = G.t and type vertex = G.V.t
 (** Classic Persistent Graphs *)
 
-module I (G : Sig.I with type V.label = int) : S with type graph = G.t
+module I (G : Sig.I with type V.label = int) :
+   S with type graph = G.t and type vertex = G.V.t
 (** Classic Imperative Graphs *)
diff --git a/src/components.ml b/src/components.ml
index 8efa43b..c727f2e 100644
--- a/src/components.ml
+++ b/src/components.ml
@@ -96,6 +96,78 @@ module Make(G: G) = struct
 
 end
 
+(** Connectivity in strongly connected directed graphs *)
+
+module Connectivity (GB: Builder.S) = struct
+
+  module MOper = Oper.Make(GB)
+  module Choose = Oper.Choose(GB.G)
+  module Dom = Dominator.Make(GB.G)
+
+  module S = Dom.S
+
+  let sstrong_articulation_points g =
+    let s = Choose.choose_vertex g in
+    let module SCC = Make (struct
+        include GB.G
+        let iter_vertex f =
+          GB.G.iter_vertex (fun v -> if not (V.equal s v) then f v)
+        let iter_succ f =
+          GB.G.iter_succ (fun v -> if not (V.equal s v) then f v)
+      end)
+    in
+    let s_is_sap = fst (SCC.scc g) > 1 in
+    let dt_s = Dom.(idom_to_dom_tree g (compute_idom g s)) in
+    let d_s = Dom.dom_tree_to_snontrivial_dom s dt_s in
+    let g_r = MOper.mirror g in
+    let dtr_s = Dom.(idom_to_dom_tree g_r (compute_idom g_r s)) in
+    let dr_s = Dom.dom_tree_to_snontrivial_dom s dtr_s in
+    let d = Dom.S.union d_s dr_s in
+    if s_is_sap then Dom.S.add s d else d
+
+  let strong_articulation_points g = S.elements (sstrong_articulation_points g)
+
+end
+
+module BiConnectivity (G: Sig.G) = struct
+
+  module Choose = Oper.Choose(G)
+  module Dom = Dominator.Make(G)
+  module RDom = Dominator.Make(
+      struct
+        type t = G.t
+        module V = G.V
+        let pred = G.succ
+        let succ = G.pred
+        let fold_vertex = G.fold_vertex
+        let iter_vertex = G.iter_vertex
+        let iter_succ = G.iter_pred
+        let nb_vertex = G.nb_vertex
+      end)
+
+  module S = Dom.S
+
+  let sstrong_articulation_points g =
+    let s = Choose.choose_vertex g in
+    let module SCC = Make (struct
+        include G
+        let iter_vertex f =
+          G.iter_vertex (fun v -> if not (V.equal s v) then f v)
+        let iter_succ f =
+          G.iter_succ (fun v -> if not (V.equal s v) then f v)
+      end)
+    in
+    let s_is_sap = fst (SCC.scc g) > 1 in
+    let dt_s = Dom.(idom_to_dom_tree g (compute_idom g s)) in
+    let d_s = Dom.dom_tree_to_snontrivial_dom s dt_s in
+    let dtr_s = RDom.(idom_to_dom_tree g (compute_idom g s)) in
+    let dr_s = Dom.dom_tree_to_snontrivial_dom s dtr_s in
+    let d = Dom.S.union d_s dr_s in
+    if s_is_sap then Dom.S.add s d else d
+
+  let strong_articulation_points g = S.elements (sstrong_articulation_points g)
+end
+
 (** Connected components (for undirected graphs) *)
 
 module type U = sig
diff --git a/src/components.mli b/src/components.mli
index ce173b5..8c98bc2 100644
--- a/src/components.mli
+++ b/src/components.mli
@@ -38,7 +38,7 @@ module Make (G: G) : sig
       number. In particular, [f u = f v] if and only if [u] and
       [v] are in the same component. Another property of the
       numbering is that components are numbered in a topological
-      order: if there is an arc from [u] to [v], then [f u >= f u]
+      order: if there is an arc from [u] to [v], then [f u >= f v]
 
       Not tail-recursive.
       Complexity: O(V+E)
@@ -56,6 +56,55 @@ module Make (G: G) : sig
 
 end
 
+(** Connectivity in strongly connected directed graphs *)
+
+module Connectivity (GB: Builder.S) : sig
+  module S : Set.S with type elt = GB.G.vertex
+
+  val strong_articulation_points : GB.G.t -> GB.G.vertex list
+  (** Computes the strong articulation points of the given
+      strongly connected directed graph. The result is undefined if the
+      input graph is not directed and strongly connected.
+
+      A strong articulation point is a vertex that when removed from the
+      original graph disconnects that graph into two or more components.
+
+      The implementation involves constructing the mirror image of the
+      graph; for bidirectional graphs prefer {!module:BiConnectivity}.
+
+      Implements the algorithm from Italiano, Laura, and Santaroni,
+      TCS 447 (2012), "Finding strong bridges and strong articulation points
+      in linear time".
+      Complexity: O(V + E) *)
+
+  val sstrong_articulation_points : GB.G.t -> S.t
+  (** As for [strong_articulation_points] but returns a set. *)
+end
+
+module BiConnectivity (G: Sig.G) : sig
+
+  module S : Set.S with type elt = G.vertex
+
+  val strong_articulation_points : G.t -> G.vertex list
+  (** Computes the strong articulation points of the given
+      strongly connected directed  graph. The result is undefined if the
+      input graph is not directed and strongly connected.
+
+      A strong articulation point is a vertex that when removed from the
+      original graph disconnects that graph into two or more components.
+
+      The implementation traverses the graph by iterating over predecessors;
+      for unidirectional graphs prefer {!module:Connectivity}.
+
+      Implements the algorithm from Italiano, Laura, and Santaroni,
+      TCS 447 (2012), "Finding strong bridges and strong articulation points
+      in linear time".
+      Complexity: O(V + E) *)
+
+  val sstrong_articulation_points : G.t -> S.t
+  (** As for [strong_articulation_points] but returns a set. *)
+end
+
 (** Connected components (for undirected graphs).
     The implementation uses union-find. Time complexity is (quasi) O(V+E).
     Space complexity is O(V). *)
diff --git a/src/cycles.ml b/src/cycles.ml
new file mode 100644
index 0000000..ea89305
--- /dev/null
+++ b/src/cycles.ml
@@ -0,0 +1,319 @@
+
+type weight =
+  | Normal of int
+  | Obligatory of int
+
+module Fashwo
+ (GB : sig
+         include Builder.S
+         val weight : G.edge -> weight
+       end)
+=
+struct
+  module G = GB.G
+
+  exception Stuck of G.vertex list
+
+  module IM = Map.Make (Int)
+  module VM = Map.Make (G.V)
+  module VS = Set.Make (G.V)
+
+  (* The algorithm of Eades, Lin, and Smyth (ELS 1993) works by "scheduling"
+     vertexes onto two lists called s1 and s2. At each iteration a vertex is
+     chosen, scheduled, and removed from the graph. Arcs from a newly scheduled
+     node toward nodes already in s1 are classified as "leftward"; they are
+     included in the generated feedback arc set. "Rightward" arcs, to vertexes
+     in s2 or that have not yet been scheduled, are not included in the
+     feedback arc set. The algorithm tries to maximize the number of rightward
+     arcs and thereby minimize the number of leftward ones. Source vertexes,
+     those with no incoming arcs in the current graph (i.e., because all its
+     predecssors have already been scheduled), are appended directly onto s1
+     and do not induce any feedback arcs. Sink vertexes are consed directly
+     onto s2 and do not induce any feedback arcs. Otherwise, the algorithm
+     chooses a vertex to maximize the difference between the number of
+     outgoing arcs and the number of incoming ones: the (remaining) incoming
+     arcs must be included in the feedback arc set. The difference between the
+     number of rightward arcs (no cost) and the number of leftward arcs
+     (feedback arcs) is called "delta". The algorithm is implemented
+     efficiently by using a data structure to group unscheduled vertexes
+     according to their delta value. When more than one vertex has the maximum
+     delta value, the original algorithm makes an arbitrary choice. The
+     algorithm of Eades and Lin (EL 1995) makes the choice using a heuristic
+     that maximizes the difference between incoming arcs and outgoing ones in
+     the vertexes that remain at the end of the iteration as such vertexes are
+     the most "unbalanced" and thus less likely to contribute to the feedback
+     arc set in future iterations. The EL 1995 algorithm includes a further
+     refinement to ignore chains of vertexes when looking for unbalanced ones,
+     since such chains do not contribute feedback arcs.
+
+     Since we just want to produce a list of feedback arcs, we don't bother
+     tracking order in s1, and we only track s2 to properly handle the
+     preprocessing optimization that removes two cycles. We maintain lists of
+     source and sink vertexes (scheduled but not yet removed from the graph)
+     and a map from delta values to sets of vertexes. As the delta value map
+     caches the state of the graph, it must be updated when the a vertex is
+     scheduled and removed from the graph. Additionally, we remember which two
+     cycles were removed during preprocessing and ensure that one of their
+     arcs is included in the feedback arc set, depending on whichever of the
+     two interlinked vertexes is scheduled first. *)
+
+  type t = {
+    s1         : VS.t;             (* vertexes placed "at left" *)
+    s2         : VS.t;             (* vertexes placed "at right";
+                                      only needed to optimize for two_cycles *)
+    sources    : VS.t;             (* vertexes with no incoming arcs *)
+    sinks      : VS.t;             (* vertexes with no outgoing arcs *)
+    delta_bins : VS.t IM.t;        (* group vertexes by delta value *)
+    vertex_bin : int VM.t;         (* map each vertex to its bin *)
+    two_cycles : G.edge list VM.t; (* edges for 2-cycles *)
+    fas        : G.edge list;      (* current feedback arc set *)
+  }
+
+  let empty = {
+      s1 = VS.empty;
+      s2 = VS.empty;
+      sources = VS.empty;
+      sinks = VS.empty;
+      delta_bins = IM.empty;
+      vertex_bin = VM.empty;
+      two_cycles = VM.empty;
+      fas = [];
+    }
+
+  let add_to_bin delta v ({ delta_bins; vertex_bin; _ } as st) =
+    { st with delta_bins =
+                IM.update delta (function None -> Some (VS.singleton v)
+                                        | Some vs -> Some (VS.add v vs))
+                  delta_bins;
+                vertex_bin = VM.add v delta vertex_bin }
+
+  let remove_from_bin v ({ delta_bins; vertex_bin; _ } as st) =
+    match VM.find_opt v vertex_bin with
+    | None -> st
+    | Some delta ->
+        { st with delta_bins =
+                    IM.update delta (function None -> None
+                                            | Some vs -> Some (VS.remove v vs))
+                      delta_bins;
+                  vertex_bin = VM.remove v vertex_bin }
+
+  (* Calculate the sums of incoming and outgoing edge weights, ignoring
+     obligatory arcs; they must be respected so their weight is irrelevant. *)
+  let weights g v =
+    let add_pweight e (s, b) =
+      match GB.weight e with Obligatory _ -> (s, true) | Normal w -> (s + w, b)
+    in
+    let add_sweight e s =
+      match GB.weight e with Obligatory w -> s + w | Normal w -> s + w
+    in
+    let inw, blocked = G.fold_pred_e add_pweight g v (0, false) in
+    let outw = G.fold_succ_e add_sweight g v 0 in
+    blocked, inw, outw
+
+  let add_vertex g v delta ({ sources; sinks; _ } as st) =
+    let ind, outd = G.in_degree g v, G.out_degree g v in
+    if ind = 0 then { st with sources = VS.add v sources }
+    else if outd = 0 then { st with sinks = VS.add v sinks }
+    else add_to_bin delta v st
+
+  (* Initialize the state for a given vertex. *)
+  let init_vertex g v st =
+    let blocked, inw, outw = weights g v in
+    if blocked then st else add_vertex g v (outw - inw) st
+
+  let init g = G.fold_vertex (init_vertex g) g empty
+
+  (* Move v from the bin for delta to sources, sinks, or another bin. *)
+  let shift_bins g v delta' st0 = add_vertex g v delta' (remove_from_bin v st0)
+
+  (* Before removing v from the graph, update the state of its sucessors. *)
+  let update_removed_succ g' e st =
+    let v = G.E.dst e in
+    let still_blocked, inw', outw' = weights g' v in
+    if still_blocked then st else shift_bins g' v (outw' - inw') st
+
+  (* Before removing v from the graph, update the state of its predecessors. *)
+  let update_removed_pred g' e ({ sinks; _ } as st) =
+    let v = G.E.src e in
+    let blocked, inw', outw' = weights g' v in
+    match GB.weight e with
+    | Obligatory _ ->
+        if blocked || outw' > 0 then st
+        else (* not blocked && outw' = 0 *)
+        { (remove_from_bin v st) with sinks = VS.add v sinks }
+    | Normal _ ->
+        if blocked then st else shift_bins g' v (outw' - inw') st
+
+  (* Remove a vertex from the graph and update the data structures for its
+     succesors and predecessors. *)
+  let remove_vertex g v st =
+    let g' = GB.remove_vertex g v in
+    (g', G.fold_succ_e (update_removed_succ g') g v st
+         |> G.fold_pred_e (update_removed_pred g') g v)
+
+  (* The original article proposes preprocessing the graph to condense long
+     chains of vertexes. This works together with the heuristic for generating
+     unbalanced vertexes, since the intermediate nodes on the chain do not
+     contribute any leftward arcs (when the last vertex is removed, they
+     become a sequence of sinks). Using such a preprocessing step with
+     weighted edges risks removing good feedback arcs, i.e., those with a big
+     difference between outgoing and incoming weights. That is why here we
+     use on-the-fly condensation, even if there is a risk of recomputing the
+     same result several times. *)
+  let rec condense w g v =
+    if G.out_degree g v = 1 then
+      match G.pred g v with
+      | [u] when not (G.V.equal u w) -> condense w g u
+      | _ -> v
+    else v
+
+  (* Find the vertex v that has the most "unbalanced" predecessor u. Most
+     unbalanced means the biggest difference between the input weights and
+     output weights. Skip any vertex with an incoming obligatory arc. *)
+  let takemax g v imax =
+    let check_edge e max = (* check u -> v *)
+      let u_blocked, u_inw, u_outw =
+        weights g (condense (G.E.dst e) g (G.E.src e)) in
+      let u_w = u_inw - u_outw in
+      match max with
+      | Some (None, _)
+      | None -> Some ((if u_blocked then None else Some u_w), v)
+      | Some (Some x_w, _) when u_w > x_w -> Some (Some u_w, v)
+      | _ -> max
+    in
+    G.fold_pred_e check_edge g v imax
+
+  (* Look for the vertex with the highest delta value that is not the target
+     of an obligatory arc. Use the "unbalanced" heuristic impllemented in
+     [takemax] to discriminate between competing possibilities. If a vertex
+     is found, remove it from the returned delta bins. *)
+  let max_from_deltas g ({ delta_bins; _ } as st) =
+    let rec f = function
+      | Seq.Nil -> None
+      | Seq.Cons ((_, dbin), tl) ->
+          (match VS.fold (takemax g) dbin None with
+           | None -> f (tl ())
+           | Some (_, v) -> Some (v, remove_from_bin v st))
+    in
+    f (IM.to_rev_seq delta_bins ())
+
+  (* Include any leftward arcs due to the two-cycles that were removed by
+     preprocessing. *)
+  let add_from_two_cycles s1 s2 two_cycles v fas =
+    let bf es b = if G.V.equal (G.E.dst b) v then b::es else es in
+    let f es e =
+      let w = G.E.dst e in
+      if VS.mem w s1 then e::es
+      else if VS.mem w s2 then
+        (* the two-cycle partner has already been scheduled as sink, so
+           the feedback edges come from it. *)
+        match VM.find_opt w two_cycles with
+        | None -> es
+        | Some bs -> List.fold_left bf es bs
+      else es in
+    match VM.find_opt v two_cycles with
+    | None -> fas
+    | Some es -> List.fold_left f fas es
+
+  (* Shift a given vertex onto s1, and add any leftward arcs to the feedback
+     arc set. *)
+  let schedule_vertex g (v, ({ s1; s2; fas; two_cycles; _ } as st)) =
+    let add_to_fas e es = if VS.mem (G.E.src e) s1 then es else e::es in
+    (v, { st with s1 = VS.add v s1;
+                  fas = G.fold_pred_e add_to_fas g v fas
+                          |> add_from_two_cycles s1 s2 two_cycles v })
+
+  (* Take the next available vertex from, in order, sources, sinks, or the
+     highset possible delta bin. *)
+  let choose_vertex g ({ s1; s2; sources; sinks; two_cycles; fas; _ } as st0) =
+    match VS.choose_opt sources with
+    | Some v ->
+        Some (v, { st0 with sources = VS.remove v sources;
+                            sinks = VS.remove v sinks;
+                            s1 = VS.add v s1;
+                            fas = add_from_two_cycles s1 s2 two_cycles v fas })
+    | None ->
+        (match VS.choose_opt sinks with
+         | Some v ->
+             Some (v, { st0 with sinks = VS.remove v sinks;
+                                 s2 = VS.add v s2;
+                                 fas = add_from_two_cycles s1 s2 two_cycles v fas })
+         | None -> Option.map (schedule_vertex g) (max_from_deltas g st0))
+
+  let add_two_cycle_edge two_cycles e =
+    VM.update (G.E.src e) (function None -> Some [e]
+                                  | Some es -> Some (e :: es)) two_cycles
+
+  let same_weight w e =
+    match GB.weight e with
+    | Obligatory _ -> false
+    | Normal w' -> w' = w
+
+  (* For every pair of distinct vertexes A and B linked to each other by
+     edges A -ab-> B and B -ba-> A with the same weight, update the mapping
+     by linking A to ab, and B to ba, and remove the edges from the graph.
+     When A is scheduled, if B is already in s1 then the edge ab is a
+     feedback arc, and similarly for B and ba. The principle is that there
+     will be a feedback arc regardless of whether A is "scheduled" before B or
+     vice versa, therefore such cycles should not constrain vertex choices. *)
+  let remove_two_cycles g0 =
+    let f e ((g, cycles) as unchanged) =
+      match GB.weight e with
+      | Obligatory _ -> unchanged
+      | Normal w ->
+          if List.length (G.find_all_edges g0 (G.E.src e) (G.E.dst e)) > 1
+          (* invalid for graphs like: { A -1-> B, A -2-> B, B -3-> A *)
+          then raise Exit
+          else
+            let back_edges =
+              G.find_all_edges g0 (G.E.dst e) (G.E.src e)
+              |> List.filter (same_weight w)
+            in
+            if back_edges = [] then unchanged
+            else (GB.remove_edge_e g e,
+                  List.fold_left add_two_cycle_edge cycles back_edges)
+    in
+    try
+      G.fold_edges_e f g0 (g0, VM.empty)
+    with Exit -> (g0, VM.empty)
+
+  (* All self loops must be broken, so just add them straight into the
+     feedback arc set. *)
+  let remove_self_loops g0 =
+    let f v (g, fas) =
+      let self_loops = G.find_all_edges g0 v v in
+      (List.fold_left GB.remove_edge_e g self_loops,
+       List.rev_append self_loops fas)
+    in
+    G.fold_vertex f g0 (g0, [])
+
+  (* Remove any arcs between strongly connected components. There can be no
+     cycles between distinct sccs by definition. *)
+  module C = Components.Make(G)
+  module Emap = Gmap.Edge(G)(struct include GB.G include GB end)
+
+  let disconnect_sccs g =
+    let nsccs, fscc = C.scc g in
+    let in_same_scc e =
+      if fscc (G.E.src e) = fscc (G.E.dst e) then Some e else None
+    in
+    if nsccs < 2 then g
+    else Emap.filter_map in_same_scc g
+
+  let feedback_arc_set g0 =
+    let rec loop (g, st) =
+      match choose_vertex g st with
+      | Some (v, st') when G.mem_vertex g v -> loop (remove_vertex g v st')
+      | Some (_, st') -> loop (g, st')
+      | None ->
+          let remaining = IM.fold (Fun.const VS.union) st.delta_bins VS.empty in
+          if VS.is_empty remaining then st.fas
+          else raise (Stuck (VS.elements remaining))
+    in
+    let g1 = disconnect_sccs g0 in
+    let g2, fas = remove_self_loops g1 in
+    let g3, two_cycles = remove_two_cycles g2 in
+    loop (g3, { (init g3) with fas; two_cycles })
+
+end
+
diff --git a/src/cycles.mli b/src/cycles.mli
new file mode 100644
index 0000000..8a47e49
--- /dev/null
+++ b/src/cycles.mli
@@ -0,0 +1,71 @@
+(**************************************************************************)
+(*                                                                        *)
+(*  Ocamlgraph: a generic graph library for OCaml                         *)
+(*  Copyright (C) 2004-2022                                               *)
+(*  Sylvain Conchon, Jean-Christophe Filliatre and Julien Signoles        *)
+(*                                                                        *)
+(*  This software is free software; you can redistribute it and/or        *)
+(*  modify it under the terms of the GNU Library General Public           *)
+(*  License version 2.1, with the special exception on linking            *)
+(*  described in file LICENSE.                                            *)
+(*                                                                        *)
+(*  This software is distributed in the hope that it will be useful,      *)
+(*  but WITHOUT ANY WARRANTY; without even the implied warranty of        *)
+(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                  *)
+(*                                                                        *)
+(**************************************************************************)
+
+(** Algorithms related to cycles in directed graphs. *)
+
+type weight =
+  | Normal of int
+    (** Weighted arc that can be included in the feedback set. The
+        weight must be zero (not normally a good choice) or positive
+        (1 may be a good choice). *)
+  | Obligatory of int
+    (** Obligatory arc that cannot be returned in the feedback set.
+        Set the weight to zero to completely ignore obligatory arcs
+        when choosing which vertex to schedule. Set it to a positive
+        value (1 may be a good choice) to adjust the preference for
+        choosing vertexes that may "unblock" other vertexes by
+        removing their incoming obligatory arcs. *)
+
+(** Adaptation of the FASH algorithm of Eades and Lin (1995) to handle
+    edge weights and obligatory arcs. The algorithm tries to minimize the
+    total weight of the returned feedback arc set. Obligatory arcs are
+    respected and never returned in the feedback arc set, an exception is
+    raised if the obligatory arcs form a cycle. The adapted algorithm is
+    hereby called FASHWO: “feedback arc set heuristic + weights and
+    obligations”.
+
+    For a graph G and any one of its feedback arc sets F, the graph G - F is
+    obviously acyclic. If F is minimal, i.e., adding any of its edges to G - F
+    would introduce a cycle, then reversing, rather than removing, the
+    feedback arcs also gives an ayclic graph, [G - F + F^R]. In fact, Eades
+    and Lin define the feedback arc set as "a set of arcs whose reversal makes
+    G acyclic".
+
+    @see <https://mathoverflow.net/a/234023/> David Epstein proof about reversed arcs *)
+module Fashwo
+ (GB : sig
+         include Builder.S
+
+         (** Assign weights to edges. *)
+         val weight : G.edge -> weight
+       end) :
+sig
+  (** Raised if cycles remain and all the remaining vertexes are obligatory.
+      The argument gives the list of remaining vertexes. *)
+  exception Stuck of GB.G.vertex list
+
+  (** Return a minimal set of arcs whose removal or reversal would make the
+      given graph acyclic.
+
+      By minimal, we mean that each arc in the returned list must be removed
+      or reversed, i.e., none are superfluous. Since a heuristic is used, the
+      returned list may not be a minimum feedback arc set. Finding the {i
+      minimum feedback arc set}, dually, the {i maximum acyclic subgraph} is
+      NP-hard for general graphs. *)
+  val feedback_arc_set : GB.G.t -> GB.G.edge list
+end
+
diff --git a/src/dominator.ml b/src/dominator.ml
index 7d6f9ff..da9dffa 100644
--- a/src/dominator.ml
+++ b/src/dominator.ml
@@ -71,6 +71,8 @@ module type S = sig
   val compute_dom_frontier: t -> dom_tree -> idom -> vertex -> vertex list
   val idom_to_dominators: ('a -> 'a) -> 'a -> 'a list
   val idom_to_dom: (vertex -> vertex) -> vertex -> vertex -> bool
+  val dom_tree_to_nontrivial_dom : vertex -> dom_tree -> vertex list
+  val dom_tree_to_snontrivial_dom : vertex -> dom_tree -> S.t
 end
 
 module Make(G : G) = struct
@@ -349,6 +351,30 @@ module Make(G : G) = struct
     with Not_found ->
       false
 
+  (* There is a nice description of non-trivial dominators with an example
+     in Section 2 and Figure 2 of Jaberi 2016, "On computing the
+     2-vertex-connected components of directed graphs". *)
+
+  let dom_tree_to_nontrivial_dom v dt =
+    let rec f rs = function
+      | [] -> rs
+      | x::xs ->
+          (match dt x with
+           | [] -> f rs xs
+           | ys -> f (x::rs) (List.rev_append ys xs))
+    in
+    f [] (dt v)
+
+  let dom_tree_to_snontrivial_dom v dt =
+    let rec f rs = function
+      | [] -> rs
+      | x::xs ->
+          (match dt x with
+           | [] -> f rs xs
+           | ys -> f (S.add x rs) (List.rev_append ys xs))
+    in
+    f S.empty (dt v)
+
 end
 
 module Make_graph(G: I) = struct
diff --git a/src/dominator.mli b/src/dominator.mli
index e3c99f8..3d91d90 100644
--- a/src/dominator.mli
+++ b/src/dominator.mli
@@ -122,6 +122,16 @@ module type S = sig
 
   val idom_to_dom: (vertex -> vertex) -> vertex -> vertex -> bool
 
+  (** Returns a list of the non-trivial dominators of a flowgraph G(v) given
+      the start vertex [v] and the corresponding dominator tree. E.g.,
+      [dom_tree_to_nontrivial_dom v (idom_to_dom_tree g (compute_idom g v))].
+      A vertex u is a non-trivial dominator of G(v) if it dominates some
+      vertex w other than v and u. *)
+  val dom_tree_to_nontrivial_dom : vertex -> dom_tree -> vertex list
+
+  (** As for [dom_tree_to_nontrivial_dom] but returns a set. *)
+  val dom_tree_to_snontrivial_dom : vertex -> dom_tree -> S.t
+
 end
 
 module Make(G : G) : S with type t = G.t and type vertex = G.V.t
diff --git a/src/eulerian.ml b/src/eulerian.ml
new file mode 100644
index 0000000..96c482c
--- /dev/null
+++ b/src/eulerian.ml
@@ -0,0 +1,269 @@
+(**************************************************************************)
+(*                                                                        *)
+(*  Ocamlgraph: a generic graph library for OCaml                         *)
+(*  Copyright (C) 2004-2010                                               *)
+(*  Sylvain Conchon, Jean-Christophe Filliatre and Julien Signoles        *)
+(*                                                                        *)
+(*  This software is free software; you can redistribute it and/or        *)
+(*  modify it under the terms of the GNU Library General Public           *)
+(*  License version 2.1, with the special exception on linking            *)
+(*  described in file LICENSE.                                            *)
+(*                                                                        *)
+(*  This software is distributed in the hope that it will be useful,      *)
+(*  but WITHOUT ANY WARRANTY; without even the implied warranty of        *)
+(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                  *)
+(*                                                                        *)
+(**************************************************************************)
+
+module type G = sig
+  type t
+  val is_directed : bool
+  module V : Sig.COMPARABLE
+  module E : Sig.EDGE with type vertex = V.t
+  val iter_edges_e : (E.t -> unit) -> t -> unit
+end
+
+(** The following implements Hierholzer's algorithm.
+
+    It is sketched as follows:
+
+    1. make a round trip from a random vertex, by following random
+       edges until we get back to the starting point (it will, as we
+       first check that all vertices have even degrees).
+
+    2. if any vertex along this cycle still has outgoing edges, pick one
+       and make another round trip from it, and then join the two cycles
+       into a single one. Repeat step 2 until all edges are exhausted.
+
+    The implementation makes use of the following:
+
+    - A table, called `out` in the following, that maps each vertex to
+      outgoing edges not yet used in the Eulerian path.
+
+    - In order to achieve optimal complexity, paths are built as
+      doubly-linked lists, so that we can merge two cycles with a common
+      vertex in constant time. This is type `dll` below.
+*)
+
+module Make(G: G) = struct
+
+  open G
+
+  let rev e =
+    E.create (E.dst e) (E.label e) (E.src e)
+
+  module H = Hashtbl.Make(V)
+
+  type out = E.t H.t H.t
+
+  let add_out_edge out x y e =
+    let s = try H.find out x
+            with Not_found -> let s = H.create 4 in H.add out x s; s in
+    H.add s y e
+
+  (** compute the table of outgoing edges *)
+  let setup g : int * out =
+    let nbe = ref 0 in
+    let out = H.create 16 in
+    let add e =
+      incr nbe;
+      let x = E.src e and y = E.dst e in
+      add_out_edge out x y e;
+      if not is_directed && not (V.equal x y) then
+        add_out_edge out y x (rev e) in
+    iter_edges_e add g;
+    !nbe, out
+
+  exception Found of V.t
+  let any h =
+    try H.iter (fun v _ -> raise (Found v)) h; assert false
+    with Found v -> v, H.find h v
+
+  type dll = { mutable prev: dll; edge: E.t; mutable next: dll }
+
+  let remove_edge out e =
+    let remove h x y =
+      let s = H.find h x in
+      assert (H.mem s y);
+      H.remove s y;
+      if H.length s = 0 then H.remove h x in
+    let v = E.src e and w = E.dst e in
+    remove out v w
+
+  let self e = V.equal (E.src e) (E.dst e)
+
+  let remove_edge edges e =
+    remove_edge edges e;
+    if not is_directed && not (self e) then remove_edge edges (rev e)
+
+  let any_out_edge out v =
+    assert (H.mem out v);
+    let s = H.find out v in
+    assert (H.length s > 0);
+    let _, e = any s in
+    remove_edge out e;
+    e
+
+  (** build an arbitrary cycle from vertex [start] *)
+  let round_trip edges start =
+    let e = any_out_edge edges start in
+    let rec path = { prev = path; edge = e; next = path } in
+    let rec tour e =
+      let v = E.dst e.edge in
+      if V.equal v start then (
+        path.prev <- e;
+        path
+      ) else (
+        let e' = { prev = e; edge = any_out_edge edges v; next = path } in
+        e.next <- e';
+        tour e'
+      ) in
+    tour path
+
+  let connect e e' =
+    e.next <- e';
+    e'.prev <- e
+
+  (** build an Eulerian cycle from vertex [start] *)
+  let eulerian_cycle out start =
+    let todos = H.create 8 in (* vertex on cycle with out edges -> cycle edge *)
+    let todo e =
+      let v = E.src e.edge in
+      if H.mem out v then H.replace todos v e else H.remove todos v in
+    let rec update start e =
+      todo e;
+      if not (V.equal (E.dst e.edge) start) then update start e.next in
+    let path = round_trip out start in
+    update start path;
+    while H.length todos > 0 do
+      let v, e = any todos in
+      H.remove todos v;
+      assert (H.mem out v);
+      let e' = round_trip out v in
+      update v e';
+      let p = e.prev in
+      assert (p.next == e);
+      let p' = e'.prev in
+      assert (p'.next == e');
+      connect p e';
+      connect p' e;
+    done;
+    path
+
+  let list_of path =
+    let rec convert acc e =
+      if e == path then List.rev acc else convert (e.edge :: acc) e.next in
+    convert [path.edge] path.next
+
+  let mem_edge out x y =
+    try H.mem (H.find out x) y with Not_found -> false
+
+  let out_degree out x =
+    try H.length (H.find out x) with Not_found -> 0
+
+  let undirected g =
+    let nbe, out = setup g in
+    let odds = H.create 2 in
+    let check v s =
+      let d = H.length s in
+      let d = if H.mem s v then d - 1 else d in
+      if d mod 2 = 1 then H.add odds v () in
+    H.iter check out;
+    let n = H.length odds in
+    if n <> 0 && n <> 2 then invalid_arg "Eulerian.path (bad degrees)";
+    let cycle = n = 0 in
+    let path =
+      if cycle then
+        if nbe = 0 then []
+        else let v, _ = any out in list_of (eulerian_cycle out v)
+      else (
+        (* we have two vertices x and y with odd degrees *)
+        let x, _ = any odds in
+        H.remove odds x;
+        let y, _ = any odds in
+
+        if mem_edge out x y then (
+          (* there is an edge x--y => it connects 1 or 2 Eulerian cycles *)
+          let xy = H.find (H.find out x) y in
+          remove_edge out xy;
+          match out_degree out x, out_degree out y with
+          | 0, 0 -> [xy]
+          | _, 0 -> rev xy :: list_of (eulerian_cycle out x)
+          | 0, _ -> xy :: list_of (eulerian_cycle out y)
+          | _ ->
+              let py = eulerian_cycle out y in
+              (* caveat: the cycle from y may exhaust edges from x *)
+              if out_degree out x = 0 then xy :: list_of py
+              else list_of (eulerian_cycle out x) @ xy :: list_of py
+                (* a bit of a pity to use list concatenation here,
+                   but this does not change the complexity *)
+        ) else (
+          (* no edge x--y => add one, build a cycle, then remove it *)
+          let dummy = E.label (snd (any (H.find out x))) in
+          let xy = E.create x dummy y in
+          H.add (H.find out x) y xy;
+          H.add (H.find out y) x (rev xy);
+          let p = eulerian_cycle out x in
+          let rec find e = (* lookup for x--y, to break the cycle there *)
+            let v = E.src e.edge and w = E.dst e.edge in
+            if V.equal v x && V.equal w y ||
+               V.equal v y && V.equal w x then e else find e.next in
+          let start = find p in
+          List.tl (list_of start)
+        )
+      )
+    in
+    (* check that all edges have been consumed *)
+    if H.length out > 0 then invalid_arg "Eulerian.path (not connected)";
+    path, cycle
+
+  let directed g =
+    let delta = H.create 16 in (* out - in *)
+    let add v d =
+      H.replace delta v (d + try H.find delta v with Not_found -> 0) in
+    let add e =
+      add (E.src e) 1; add (E.dst e) (-1) in
+    iter_edges_e add g;
+    let start = ref None and finish = ref None in
+    let check v = function
+      | 1 when !start = None -> start := Some v
+      | -1 when !finish = None -> finish := Some v
+      | 0 -> ()
+      | _ -> invalid_arg "Eulerian.path (bad degrees)" in
+    H.iter check delta;
+    let nbe, out = setup g in
+    let path, cycle = match !start, !finish with
+      | None, None when nbe = 0 ->
+          [], true
+      | None, None ->
+          let v, _ = any out in list_of (eulerian_cycle out v), true
+      | Some s, Some f ->
+          (* add one edge f->s, build a cycle, then remove it
+             note: there may be already an edge f->s
+                   if so, we are adding *a second one* and we are careful
+                   about removing this one, not the other *)
+          let dummy = E.label (snd (any (H.find out s))) in
+          let fs = E.create f dummy s in
+          add_out_edge out f s fs;
+          let p = eulerian_cycle out s in
+          let rec find e = (* lookup for f->s, to break the cycle there *)
+            if e.edge == fs then e else find e.next in
+          let start = find p in
+          List.tl (list_of start), false
+      | Some _, None
+      | None, Some _ ->
+          assert false (* since the sum of all deltas is zero *)
+    in
+    (* check that all edges have been consumed *)
+    if H.length out > 0 then invalid_arg "Eulerian.path (not connected)";
+    path, cycle
+
+  let path =
+    if is_directed then directed else undirected
+
+  let cycle g =
+    let p, c = path g in
+    if not c then invalid_arg "Eulerian.cycle";
+    p
+
+end
diff --git a/src/eulerian.mli b/src/eulerian.mli
new file mode 100644
index 0000000..640b923
--- /dev/null
+++ b/src/eulerian.mli
@@ -0,0 +1,46 @@
+(**************************************************************************)
+(*                                                                        *)
+(*  Ocamlgraph: a generic graph library for OCaml                         *)
+(*  Copyright (C) 2004-2010                                               *)
+(*  Sylvain Conchon, Jean-Christophe Filliatre and Julien Signoles        *)
+(*                                                                        *)
+(*  This software is free software; you can redistribute it and/or        *)
+(*  modify it under the terms of the GNU Library General Public           *)
+(*  License version 2.1, with the special exception on linking            *)
+(*  described in file LICENSE.                                            *)
+(*                                                                        *)
+(*  This software is distributed in the hope that it will be useful,      *)
+(*  but WITHOUT ANY WARRANTY; without even the implied warranty of        *)
+(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                  *)
+(*                                                                        *)
+(**************************************************************************)
+
+(** Eulerian path
+
+    This module implements Hierholzer's algorithm, with O(E) complexity
+    where E is the number of edges.
+
+    Limitations:
+    - multigraphs are not supported
+ *)
+
+module type G = sig
+  type t
+  val is_directed : bool
+  module V : Sig.COMPARABLE
+  module E : Sig.EDGE with type vertex = V.t
+  val iter_edges_e : (E.t -> unit) -> t -> unit
+end
+
+module Make(G: G) : sig
+
+  val path: G.t -> G.E.t list * bool
+  (** [path g] returns an Eulerian path of [g]. The Boolean indicates
+      whether the path is a cycle. Raises [Invalid_argument] if there
+      is no Eulerian path. *)
+
+  val cycle: G.t -> G.E.t list
+  (** [cycle g] returns an Eulerian cycle of [g].
+      Raises [Invalid_argument] if there is no Eulerian cycle. *)
+
+end
diff --git a/src/graph.ml b/src/graph.ml
index 08bc93a..5f374dc 100644
--- a/src/graph.ml
+++ b/src/graph.ml
@@ -12,6 +12,7 @@ module Rand = Rand
 module Oper = Oper
 module Components = Components
 module Path = Path
+module Cycles = Cycles
 module Nonnegative = Nonnegative
 module Traverse = Traverse
 module Coloring = Coloring
diff --git a/src/graphviz.ml b/src/graphviz.ml
index 0d05da7..75ce965 100644
--- a/src/graphviz.ml
+++ b/src/graphviz.ml
@@ -576,7 +576,7 @@ struct
     fprintf ppf "%t@ " print_nodes;
     fprintf ppf "%t@ " print_subgraphs;
     fprintf ppf "%t@ " print_edges;
-    fprintf ppf "@]}@]"
+    fprintf ppf "@]}@;@]"
 
   (** [output_graph oc graph] pretty prints the graph [graph] in the dot
       language on the channel [oc]. *)
diff --git a/src/imperative.ml b/src/imperative.ml
index 0044da3..96fbd94 100644
--- a/src/imperative.ml
+++ b/src/imperative.ml
@@ -477,15 +477,13 @@ module Matrix = struct
     (* map iterator on vertex *)
     let map_vertex f g =
       let n = nb_vertex g in
+      let f i = (* ensures f is applied exactly once for each vertex *)
+        let fi = f i in
+        if fi < 0 || fi >= n then invalid_arg "[ocamlgraph] map_vertex";
+        fi in
+      let v = Array.init n f in
       let g' = make n in
-      iter_edges
-        (fun i j ->
-           let fi = f i in
-           let fj = f j in
-           if fi < 0 || fi >= n || fj < 0 || fj >= n then
-             invalid_arg "[ocamlgraph] map_vertex";
-           Bitv.unsafe_set g'.(fi) fj true)
-        g;
+      iter_edges (fun i j -> Bitv.unsafe_set g'.(v.(i)) v.(j) true) g;
       g'
 
     (* labeled edges going from/to a vertex *)
diff --git a/src/pack.ml b/src/pack.ml
index 4b8c14d..0436ef4 100644
--- a/src/pack.ml
+++ b/src/pack.ml
@@ -93,6 +93,10 @@ struct
     let iter_stable = S.iter
   end
 
+  module Eulerian = struct
+    include Eulerian.Make(G)
+  end
+
   module Int = struct
     type t = int
     let compare : t -> t -> int = Stdlib.compare
diff --git a/src/path.ml b/src/path.ml
index fb5e0cd..aebf4c4 100644
--- a/src/path.ml
+++ b/src/path.ml
@@ -327,6 +327,9 @@ struct
       (* the path is not in cache; we check it with a BFS *)
       let visited = HV.create 97 in
       let q = Queue.create () in
+      (* [visited] contains exactly the vertices that have been added to [q] *)
+      let push v =
+        if not (HV.mem visited v) then (HV.add visited v (); Queue.add v q) in
       let rec loop () =
         if Queue.is_empty q then begin
           HVV.add pc.cache (v1, v2) false;
@@ -337,15 +340,12 @@ struct
           if G.V.compare v v2 = 0 then
             true
           else begin
-            if not (HV.mem visited v) then begin
-              HV.add visited v ();
-              G.iter_succ (fun v' -> Queue.add v' q) pc.graph v
-            end;
+            G.iter_succ push pc.graph v;
             loop ()
           end
         end
       in
-      Queue.add v1 q;
+      push v1;
       loop ()
 
 end
diff --git a/src/path.mli b/src/path.mli
index 4c1f222..480d4ef 100644
--- a/src/path.mli
+++ b/src/path.mli
@@ -142,8 +142,8 @@ sig
 
       Complexity: The path checker contains a cache of all results computed
       so far. This cache is implemented with a hash table so access in this
-      cache is usually O(1). When the result is not in the cache, Dijkstra's
-      algorithm is run to check for the path, and all intermediate results
+      cache is usually O(1). When the result is not in the cache, a BFS
+      is run to check for the path, and all intermediate results
       are cached.
 
       Note: if checks are to be done for almost all pairs of vertices, it
diff --git a/src/sig.mli b/src/sig.mli
index abe28f9..b77e8cb 100644
--- a/src/sig.mli
+++ b/src/sig.mli
@@ -189,7 +189,13 @@ module type G = sig
   (** Fold on all edges of a graph. *)
 
   val map_vertex : (vertex -> vertex) -> t -> t
-  (** Map on all vertices of a graph. *)
+  (** Map on all vertices of a graph.
+
+      The current implementation requires the supplied function to be
+      injective. Said otherwise, [map_vertex] cannot be used to contract
+      a graph by mapping several vertices to the same vertex.
+      To contract a graph, use instead [create], [add_vertex],
+      and [add_edge]. *)
 
   (** {2 Vertex iterators}
 
diff --git a/src/sig_pack.mli b/src/sig_pack.mli
index 68fd77f..143d493 100644
--- a/src/sig_pack.mli
+++ b/src/sig_pack.mli
@@ -344,6 +344,18 @@ module type S = sig
     (** [full n] builds a graph with [n] vertices and all possible edges.
         The optional argument [self] indicates if loop edges should be added
         (default value is [true]). *)
+
+    val cycle : int -> t * V.t array
+    (** [cycle n] builds a graph that is a cycle with [n] vertices.
+        Vertices are labelled with 0,1,...,n-1 and there is an edge from
+        vertex [i] to vertex [(i+1) mod n].
+        Vertices are also returned in an array for convenience. *)
+
+    val grid : n:int -> m:int -> t * V.t array array
+    (** [grid n m] builds a grid graph with [n*m] vertices, with edges
+        from vertex [(i,j)] to vertices [(i+1,j)] and [(i,j+1)] (and no
+        wrapping around). Vertex [(i,j)] is labelled with [i*m+j].
+        Vertices are also returned in a [n*m] matrix for convenience. *)
   end
 
   (** Random graphs *)
@@ -409,6 +421,16 @@ module type S = sig
     val iter_stable : (V.t -> unit) -> t -> unit
   end
 
+  (** Eulerian path *)
+  module Eulerian : sig
+    val path: t -> E.t list * bool
+      (** [path g] returns an Eulerian path of g. The Boolean indicates
+          whether the path is a cycle. Raises [Invalid_argument] if there is
+          no Eulerian path. *)
+
+    val cycle: t -> E.t list
+  end
+
   val spanningtree : t -> E.t list
   (** Kruskal algorithm *)
 
diff --git a/src/traverse.mli b/src/traverse.mli
index 122f877..975be54 100644
--- a/src/traverse.mli
+++ b/src/traverse.mli
@@ -17,7 +17,22 @@
 
 (** Graph traversal. *)
 
-(** {2 Dfs and Bfs} *)
+(** {2 Dfs and Bfs}
+
+   In the modules below, the most meaningful functions are the
+   [iter/fold_component] functions, where the starting point of the
+   traversal is user-provided.
+
+   Functions [iter/fold] to traverse the whole graph are also
+   provided, for convenience, and they proceed as follows: they run
+   the user-provided [iter/fold_vertex] functions (from input module
+   [G]) and, for each vertex not yet visited, start a new traversal
+   from this vertex. In particular, each traversal is not necessarily
+   started from a vertex without predecessors. Said otherwise, it is
+   up to you to come up with an [iter_vertex] function that will
+   identify suitable roots, e.g. vertices with no predecessors, if
+   this is really what you want.
+*)
 
 (** Minimal graph signature for {!Dfs} and {!Bfs}.
     Sub-signature of {!Sig.G}. *)
@@ -27,11 +42,13 @@ module type G = sig
   module V : Sig.COMPARABLE
   val iter_vertex : (V.t -> unit) -> t -> unit
   (** It is enough to iter over all the roots (vertices without predecessor) of
-      the graph, even if iterating over the other vertices is correct. *)
+      the graph, even if iterating over the other vertices is correct.
+      (See the comment above.) *)
 
   val fold_vertex : (V.t -> 'a -> 'a) -> t  -> 'a -> 'a
   (** It is enough to fold over all the roots (vertices without predecessor) of
-      the graph, even if folding over the other vertices is correct. *)
+      the graph, even if folding over the other vertices is correct.
+      (See the comment above.) *)
 
   val iter_succ : (V.t -> unit) -> t -> V.t -> unit
   val fold_succ : (V.t -> 'a -> 'a) -> t -> V.t -> 'a -> 'a
@@ -68,7 +85,7 @@ module Dfs(G : G) : sig
 
   val fold : (G.V.t -> 'a -> 'a) -> 'a -> G.t -> 'a
   (** The function is applied each time a node is reached for the first time,
-      before idoterating over its successors. Tail-recursive. *)
+      before iterating over its successors. Tail-recursive. *)
 
   val fold_component : (G.V.t -> 'a -> 'a) -> 'a -> G.t -> G.V.t -> 'a
   (** Idem, but limited to a single root vertex. *)
@@ -102,12 +119,20 @@ module Bfs(G : G) : sig
   (** {2 Classical big-step iterators} *)
 
   val iter : (G.V.t -> unit) -> G.t -> unit
+  (** The function is applied each time a node is reached for the first time.
+      Not tail-recursive. *)
+
   val iter_component : (G.V.t -> unit) -> G.t -> G.V.t -> unit
+  (** Idem, but limited to a single root vertex. *)
 
   (** {2 Classical folds} *)
 
   val fold : (G.V.t -> 'a -> 'a) -> 'a -> G.t -> 'a
+  (** The function is applied each time a node is reached for the first time.
+      Not tail-recursive. *)
+
   val fold_component : (G.V.t -> 'a -> 'a) -> 'a -> G.t -> G.V.t -> 'a
+  (** Idem, but limited to a single root vertex. *)
 
   (** {2 Step-by-step iterator}
       See module [Dfs] *)
diff --git a/src/util.ml b/src/util.ml
index a7d50ea..afbfb36 100644
--- a/src/util.ml
+++ b/src/util.ml
@@ -48,3 +48,9 @@ module DataV(L : sig type t end)(V : Sig.COMPARABLE) = struct
   let set_data (y, _) = (:=) y
 end
 
+module Memo(X: HASHABLE) = struct
+  module H = Hashtbl.Make(X)
+  let memo ?(size=128) f =
+    let h = H.create size in
+    fun x -> try H.find h x with Not_found -> let y = f x in H.add h x y; y
+end
diff --git a/src/util.mli b/src/util.mli
index 63ae42f..165296b 100644
--- a/src/util.mli
+++ b/src/util.mli
@@ -45,3 +45,6 @@ module DataV(L : sig type t end)(V : Sig.COMPARABLE) : sig
   val set_data : t -> data -> unit
 end
 
+module Memo(X: HASHABLE) : sig
+  val memo: ?size:int -> (X.t -> 'a) -> X.t -> 'a
+end
diff --git a/tests/dot.expected b/tests/dot.expected
index fd3164b..4fb53d9 100644
--- a/tests/dot.expected
+++ b/tests/dot.expected
@@ -14,4 +14,5 @@ digraph G {
   "task_cook" -> "task_eat" [label=<f&#36;oo>, ];
   "task_table" -> "task_eat" [label=<f&#36;oo>, ];
   
-  }========= END output graph =========
\ No newline at end of file
+  }
+========= END output graph =========
\ No newline at end of file
diff --git a/tests/dune b/tests/dune
index 12a0d45..840bf7b 100644
--- a/tests/dune
+++ b/tests/dune
@@ -8,6 +8,21 @@
  (libraries graph)
  (modules test_topsort))
 
+(test
+ (name test_check_path)
+ (libraries graph)
+ (modules test_check_path))
+
+(test
+ (name test_map_vertex)
+ (libraries graph)
+ (modules test_map_vertex))
+
+(test
+ (name test_eulerian)
+ (libraries graph)
+ (modules test_eulerian))
+
 ;; Rules for the Bellman-Ford tests
 
 (rule
@@ -103,6 +118,63 @@
  (modules test_johnson)
  (libraries graph))
 
+;; Rules for the test_nontrivial_dom test
+
+(rule
+ (with-stdout-to
+  test_nontrivial_dom.output
+  (run ./test_nontrivial_dom.exe)))
+
+(rule
+ (alias runtest)
+ (action
+  (progn
+   (diff test_nontrivial_dom.expected test_nontrivial_dom.output)
+   (echo "test_nontrivial_dom: all tests succeeded.\n"))))
+
+(executable
+ (name test_nontrivial_dom)
+ (modules test_nontrivial_dom)
+ (libraries graph))
+
+;; Rules for the test_saps test
+
+(rule
+ (with-stdout-to
+  test_saps.output
+  (run ./test_saps.exe)))
+
+(rule
+ (alias runtest)
+ (action
+  (progn
+   (diff test_saps.expected test_saps.output)
+   (echo "test_saps: all tests succeeded.\n"))))
+
+(executable
+ (name test_saps)
+ (modules test_saps)
+ (libraries graph))
+
+;; Rules for the test_cycles test
+
+(rule
+ (with-stdout-to
+  test_cycles.output
+  (run ./test_cycles.exe)))
+
+(rule
+ (alias runtest)
+ (action
+  (progn
+   (diff test_cycles.expected test_cycles.output)
+   (echo "test_cycles: all tests succeeded.\n"))))
+
+(executable
+ (name test_cycles)
+ (modules test_cycles)
+ (libraries graph))
+
 ;; Rules for the weak topological test
 
 (rule
diff --git a/tests/test_31.ml b/tests/test_31.ml
new file mode 100644
index 0000000..2c59368
--- /dev/null
+++ b/tests/test_31.ml
@@ -0,0 +1,32 @@
+
+(* Issue #31
+
+   This would fail in the current state of OCamlGraph,
+   as map_vertex requires a supplied function that is injective.
+*)
+
+module String = struct include String let hash = Hashtbl.hash end
+module G = Graph.Persistent.Digraph.ConcreteBidirectional(String)
+
+(* Make a persistent graph where:
+     A -> B
+     B -> C
+     C -> B *)
+let g = List.fold_left (fun g (x, y) ->
+            G.add_edge g x y) G.empty [("A", "B"); ("B", "C"); ("C", "B")]
+
+(* Contract the SCC *)
+module C = Graph.Components.Make(G)
+
+let contract g =
+  let names = Array.map (String.concat ",") (C.scc_array g) in
+  let (_, numbering) = C.scc g in
+  Array.fold_left (fun g x -> G.remove_edge g x x)
+    (G.map_vertex (fun x -> names.(numbering x)) g)
+    names
+
+let g' = contract g
+(* this is now "A" -> "C,B" *)
+
+let () =
+  assert (G.in_degree g' "C,B" = 1)
diff --git a/tests/test_check_path.ml b/tests/test_check_path.ml
new file mode 100644
index 0000000..b5cc4ac
--- /dev/null
+++ b/tests/test_check_path.ml
@@ -0,0 +1,52 @@
+
+(* Test file for Path.Check *)
+
+open Format
+open Graph
+open Pack.Digraph
+
+let test n edges =
+  let v = Array.init n V.create in
+  let g = create () in
+  let () = Array.iter (add_vertex g) v in
+  let build (s,t) = add_edge g v.(s) v.(t) in
+  List.iter build edges;
+  let path = PathCheck.check_path (PathCheck.create g) in
+  for i = 0 to n - 1 do
+    let seen = Array.make n false in
+    let pre v = seen.(V.label v) <- true in
+    Dfs.prefix_component pre g v.(i);
+    for j = 0 to n - 1 do
+      assert (seen.(j) = path v.(i) v.(j))
+    done
+  done
+
+let () =
+  test 3 [0,1; 1,2];
+  test 3 [];
+  (* 1-cycle *)
+  test 1 [0,0];
+  (* 2-cycle *)
+  test 2 [0,1; 1,0];
+  test 3 [0,1; 1,0];
+  (* 2-cycle with out edge *)
+  test 3 [0,1; 1,0; 1,2];
+  test 3 [2,0; 0,2; 0,1];
+  test 3 [1,2; 2,1; 2,0];
+  (* 2 loops *)
+  test 5 [1,2; 2,1; 2,0; 3,4; 4,3];
+  test 5 [1,2; 2,1; 2,0; 2,3; 3,4; 4,3];
+  (* 2-cycle with in edge *)
+  test 3 [1,2; 2,1; 0,2];
+  test 3 [1,2; 2,1; 0,1];
+  (* 2 cycles connected *)
+  test 4 [0,1; 1,0; 2,3; 3,2; 2,1];
+  test 4 [0,1; 1,0; 2,3; 3,2; 1,2];
+  test 4 [0,1; 1,0; 2,3; 3,2; 1,2; 2,1];
+  (* 3-cycle with in and out edges *)
+  test 5 [0,1; 1,2; 2,0; 3,0; 2,4];
+  (* 3 cycles in a row *)
+  test 7 [0,1; 1,0; 1,2; 2,3; 3,2; 3,4; 4,5; 5,6; 6,4];
+  (* 3 cycles with 2 cycles in a cycle *)
+  test 7 [0,1; 1,0; 1,2; 2,3; 3,2; 3,4; 4,5; 5,6; 6,4; 5,2];
+  printf "test check_path: all tests succeeded.@."
diff --git a/tests/test_components.expected b/tests/test_components.expected
index fe0501b..dd46d10 100644
--- a/tests/test_components.expected
+++ b/tests/test_components.expected
@@ -1,11 +1,11 @@
 7 components
 0 -> 0
 1 -> 1
-2 -> 1
-3 -> 2
-4 -> 3
-5 -> 3
+2 -> 2
+3 -> 3
+4 -> 1
+5 -> 1
 6 -> 4
-7 -> 5
-8 -> 6
-9 -> 0
+7 -> 1
+8 -> 5
+9 -> 6
diff --git a/tests/test_components.ml b/tests/test_components.ml
index 222d407..f052a57 100644
--- a/tests/test_components.ml
+++ b/tests/test_components.ml
@@ -23,6 +23,7 @@ module C = Components.Undirected(Pack.Graph)
 open Pack.Graph
 
 let () =
+  Random.init 42;
   let g = Rand.graph ~v:10 ~e:3 () in
   let n, f = C.components g in
   printf "%d components@." n;
diff --git a/tests/test_cycles.expected b/tests/test_cycles.expected
new file mode 100644
index 0000000..869bd71
--- /dev/null
+++ b/tests/test_cycles.expected
@@ -0,0 +1,3 @@
+cycles1 = { 3 -> 2 } (cycles to no cycles)
+cycles2 = { 7 -> 5, 3 -> 1 } (cycles to no cycles)
+cycles3 = { 4 -> 2, 3 -> 1 } (cycles to no cycles)
diff --git a/tests/test_cycles.ml b/tests/test_cycles.ml
new file mode 100644
index 0000000..5765bb2
--- /dev/null
+++ b/tests/test_cycles.ml
@@ -0,0 +1,100 @@
+
+(* Test file for Cycles module *)
+
+open Graph
+
+module Int = struct
+  type t = int
+  let compare = compare
+  let hash = Hashtbl.hash
+  let equal = (=)
+  let default = 0
+end
+
+let pp_comma p () = Format.(pp_print_char p ','; pp_print_space p ())
+let pp_edge p (s, d) = Format.fprintf p "%d -> %d" s d
+
+module GP = Persistent.Digraph.Concrete(Int)
+
+module GPDFS = Traverse.Dfs (GP)
+
+let pp_has_cycles p g =
+  if GPDFS.has_cycle g
+  then Format.pp_print_string p "cycles"
+  else Format.pp_print_string p "no cycles"
+
+module FW = Cycles.Fashwo(struct
+    include Builder.P(GP)
+    let weight _ = Cycles.Normal 1
+  end)
+
+(* Eades and Linh, "A Heuristic for the Feedback Arc Set Problem", Fig. 1 *)
+let g1 =
+  List.fold_left (fun g (s, d) -> GP.add_edge g s d) GP.empty
+    [ (1, 4);
+      (1, 3);
+      (2, 1);
+      (2, 4);
+      (3, 2);
+      (4, 3);
+    ]
+let cycles1 = FW.feedback_arc_set g1
+let g1' = List.fold_left (fun g (s, d) -> GP.remove_edge g s d) g1 cycles1
+
+let () =
+  Format.(printf "cycles1 = @[<hv 2>{ %a }@] (%a to %a)@."
+    (pp_print_list ~pp_sep:pp_comma pp_edge) cycles1
+    pp_has_cycles g1
+    pp_has_cycles g1')
+
+(* Eades and Linh, "A Heuristic for the Feedback Arc Set Problem", Fig. 5 *)
+let g2 =
+  List.fold_left (fun g (s, d) -> GP.add_edge g s d) GP.empty
+    [ (1, 2);
+      (1, 4);
+      (2, 3);
+      (2, 4);
+      (3, 1);
+      (4, 8);
+      (5, 3);
+      (5, 6);
+      (6, 7);
+      (7, 5);
+      (8, 6);
+      (8, 7);
+    ]
+let cycles2 = FW.feedback_arc_set g2
+let g2' = List.fold_left
+            (fun g (s, d) -> GP.add_edge (GP.remove_edge g s d) d s)
+            g2 cycles2
+
+let () =
+  Format.(printf "cycles2 = @[<hv 2>{ %a }@] (%a to %a)@."
+    (pp_print_list ~pp_sep:pp_comma pp_edge) cycles2
+    pp_has_cycles g2
+    pp_has_cycles g2')
+
+(* Eades and Linh, "A Heuristic for the Feedback Arc Set Problem", Fig. 6 *)
+let g3 =
+  List.fold_left (fun g (s, d) -> GP.add_edge g s d) GP.empty
+    [ (1, 2);
+      (1, 5);
+      (2, 6);
+      (3, 1);
+      (4, 2);
+      (4, 3);
+      (5, 3);
+      (5, 6);
+      (6, 4);
+    ]
+let cycles3 = FW.feedback_arc_set g3
+let g3' = List.fold_left
+            (fun g (s, d) -> GP.add_edge (GP.remove_edge g s d) d s)
+            g3 cycles3
+
+let () =
+  Format.(printf "cycles3 = @[<hv 2>{ %a }@] (%a to %a)@."
+    (pp_print_list ~pp_sep:pp_comma pp_edge) cycles3
+    pp_has_cycles g3
+    pp_has_cycles g3')
+
diff --git a/tests/test_eulerian.ml b/tests/test_eulerian.ml
new file mode 100644
index 0000000..ee6a39b
--- /dev/null
+++ b/tests/test_eulerian.ml
@@ -0,0 +1,140 @@
+
+open Format
+open Graph
+
+open Pack.Graph
+
+let print_vertex fmt v =
+  fprintf fmt "%d" (V.label v)
+let print_edge fmt e =
+  fprintf fmt "%a->%a" print_vertex (E.src e) print_vertex (E.dst e)
+let print_path fmt p =
+  List.iter (fun e -> fprintf fmt "%a " print_edge e) p
+
+module G = Pack.Graph
+
+let exists_path g =
+  try ignore (Eulerian.path g); true with Invalid_argument _ -> false
+let exists_cycle g =
+  try ignore (Eulerian.cycle g); true with Invalid_argument _ -> false
+
+let g = create ()
+let add_vertex i = let v = V.create i in add_vertex g v; v
+let path_length g = let p, _ = Eulerian.path g in List.length p
+
+let v0 = add_vertex 0
+let () = assert (exists_path g)
+let () = assert (exists_cycle g)
+
+let v1 = add_vertex 1
+let () = assert (exists_path g)
+let () = assert (exists_cycle g)
+
+let () = add_edge g v0 v1
+let () = assert (exists_path g)
+let () = assert (not (exists_cycle g))
+let () = assert (path_length g = 1)
+
+let v2 = add_vertex 2
+let () = add_edge g v1 v2
+let () = assert (exists_path g)
+let () = assert (not (exists_cycle g))
+let () = assert (path_length g = 2)
+
+let () = add_edge g v2 v0
+let () = assert (exists_path g)
+let () = assert (exists_cycle g)
+let () = assert (path_length g = 3)
+
+let () = add_edge g v0 v0
+let () = assert (exists_cycle g)
+
+let v3 = add_vertex 3
+let () = add_edge g v2 v3
+let () = assert (exists_path g)
+let () = assert (not (exists_cycle g))
+let () = assert (path_length g = 5)
+
+let v4 = add_vertex 4
+let () = add_edge g v3 v4
+let () = add_edge g v2 v4
+let () = assert (exists_cycle g)
+let () = assert (path_length g = 7)
+
+let () = remove_edge g v2 v4
+let v5 = add_vertex 5
+let () = add_edge g v4 v5
+let () = add_edge g v5 v3
+let () = assert (exists_path g)
+let () = assert (not (exists_cycle g))
+let () = assert (path_length g = 8)
+
+let () = remove_edge g v2 v3 (* not connected anymore *)
+let () = assert (not (exists_path g))
+
+let () =
+  for n = 2 to 5 do
+    let g = Classic.full ~self:false  (2*n) in
+    assert (not (exists_path g));
+    let g = Classic.full ~self:true  (2*n) in
+    assert (not (exists_path g));
+    let g = Classic.full ~self:false (2*n+1) in
+    let p, c = Eulerian.path g in
+    assert c;
+    assert (List.length p = n*(2*n+1));
+    let g = Classic.full ~self:true (2*n+1) in
+    let p, c = Eulerian.path g in
+    assert c;
+    assert (List.length p = (n+1)*(2*n+1))
+   done
+
+let () =
+  (* +---x---+  tricky one, as the edge x-y
+     |   |   |  connects two vertices on a cycle
+     +---y---+ *)
+  let g, _ = Classic.grid ~n:2 ~m:3 in
+  let p, c = Eulerian.path g in
+  assert (not c);
+  assert (List.length p = 7)
+
+open Pack.Digraph
+
+let exists_path g =
+  try ignore (Eulerian.path g); true with Invalid_argument _ -> false
+let exists_cycle g =
+  try ignore (Eulerian.cycle g); true with Invalid_argument _ -> false
+
+let () =
+  for n = 0 to 4 do
+    let g, v = Classic.cycle n in
+    let p, c = Eulerian.path g in
+    assert c;
+    assert (List.length p = n);
+    if n > 1 then (
+      remove_edge g v.(0) v.(1);
+      let p, c = Eulerian.path g in
+      assert (not c);
+      assert (List.length p = n - 1);
+    )
+  done
+
+let g, v = Classic.cycle 5
+let () = add_edge g v.(1) v.(4)
+let () = assert (not (exists_cycle g))
+let () = assert (exists_path g)
+let () = add_edge g v.(4) v.(1)
+let () = assert (exists_cycle g)
+
+(*    +------- 2 <----+
+      v               |
+    0(finish) ------> 1(start)
+      ^               |
+      +------- 3 <----+      *)
+
+let g, v = Classic.cycle 3
+let v3 = V.create 3
+let () = add_vertex g v3; add_edge g v.(1) v3; add_edge g v3 v.(0)
+let _, c = Eulerian.path g
+let () = assert (not c)
+
+
diff --git a/tests/test_map_vertex.ml b/tests/test_map_vertex.ml
new file mode 100644
index 0000000..a442494
--- /dev/null
+++ b/tests/test_map_vertex.ml
@@ -0,0 +1,80 @@
+
+(* Check that map_vertex applies the function exactly once per vertex *)
+
+open Graph
+
+let () = Random.init 1597
+
+module TestB(B: Builder.S with type G.V.label = int) = struct
+  let test n =
+    let v = Array.init n B.G.V.create in
+    let rec make g i =
+      if i = n then g else make (B.add_vertex g v.(i)) (i + 1) in
+    let g = ref (make (B.empty ()) 0) in
+    for i = 0 to n - 1 do
+      for j = 0 to n - 1 do
+        if Random.bool () then g := B.add_edge !g v.(i) v.(j)
+      done
+    done;
+    let counter = ref 0 in
+    let f x = incr counter; x in
+    let g' = B.G.map_vertex f !g in
+    assert (!counter = n);
+    assert (B.G.nb_vertex g' = n)
+
+  let () =
+    for n = 0 to 10 do test n done
+end
+module TestI(G: Sig.I with type V.label = int) = TestB(Builder.I(G))
+module TestP(G: Sig.P with type V.label = int) = TestB(Builder.P(G))
+
+module Int = struct include Int let hash x = x let default = 42 end
+
+include TestI(Pack.Digraph)
+include TestI(Pack.Graph)
+
+(* imperative, directed *)
+include TestI(Imperative.Digraph.Concrete(Int))
+include TestI(Imperative.Digraph.Abstract(Int))
+include TestI(Imperative.Digraph.ConcreteBidirectional(Int))
+include TestI(Imperative.Digraph.ConcreteLabeled(Int)(Int))
+include TestI(Imperative.Digraph.AbstractLabeled(Int)(Int))
+include TestI(Imperative.Digraph.ConcreteBidirectionalLabeled(Int)(Int))
+(* imperative, undirected *)
+include TestI(Imperative.Graph.Concrete(Int))
+include TestI(Imperative.Graph.Abstract(Int))
+include TestI(Imperative.Graph.ConcreteLabeled(Int)(Int))
+include TestI(Imperative.Graph.AbstractLabeled(Int)(Int))
+
+module TestM(G: Imperative.Matrix.S) = struct
+  let test n =
+    let g = G.make n in
+    for i = 0 to n - 1 do
+      for j = 0 to n - 1 do
+        if Random.bool () then G.add_edge g i j
+      done
+    done;
+    let counter = ref 0 in
+    let f x = incr counter; x in
+    let g' = G.map_vertex f g in
+    assert (!counter = n);
+    assert (G.nb_vertex g' = n)
+
+  let () =
+    for n = 0 to 10 do test n done
+end
+include TestM(Imperative.Matrix.Digraph)
+include TestM(Imperative.Matrix.Graph)
+
+(* persistent, directed *)
+include TestP(Persistent.Digraph.Concrete(Int))
+include TestP(Persistent.Digraph.Abstract(Int))
+include TestP(Persistent.Digraph.ConcreteBidirectional(Int))
+include TestP(Persistent.Digraph.ConcreteLabeled(Int)(Int))
+include TestP(Persistent.Digraph.AbstractLabeled(Int)(Int))
+include TestP(Persistent.Digraph.ConcreteBidirectionalLabeled(Int)(Int))
+(* persistent, undirected *)
+include TestP(Persistent.Graph.Concrete(Int))
+include TestP(Persistent.Graph.Abstract(Int))
+include TestP(Persistent.Graph.ConcreteLabeled(Int)(Int))
+include TestP(Persistent.Graph.AbstractLabeled(Int)(Int))
diff --git a/tests/test_nontrivial_dom.expected b/tests/test_nontrivial_dom.expected
new file mode 100644
index 0000000..ec64819
--- /dev/null
+++ b/tests/test_nontrivial_dom.expected
@@ -0,0 +1,2 @@
+{3, 6, 9}
+{3, 6, 9}
diff --git a/tests/test_nontrivial_dom.ml b/tests/test_nontrivial_dom.ml
new file mode 100644
index 0000000..78b70f3
--- /dev/null
+++ b/tests/test_nontrivial_dom.ml
@@ -0,0 +1,57 @@
+
+(* Test file for Dominators.dom_tree_to_(s)nontrivial_dom *)
+
+open Graph
+
+module Int = struct
+  type t = int
+  let compare = compare
+  let hash = Hashtbl.hash
+  let equal = (=)
+  let default = 0
+  end
+
+module G = Imperative.Digraph.ConcreteLabeled(Int)(Int)
+
+module Dominator = Dominator.Make (G)
+
+let g = G.create ()
+
+let () =
+  G.add_edge g 1 11;
+  G.add_edge g 1 2;
+  G.add_edge g 1 9;
+  G.add_edge g 2 3;
+  G.add_edge g 2 9;
+  G.add_edge g 2 10;
+  G.add_edge g 3 4;
+  G.add_edge g 3 5;
+  G.add_edge g 4 5;
+  G.add_edge g 5 6;
+  G.add_edge g 6 8;
+  G.add_edge g 7 6;
+  G.add_edge g 8 1;
+  G.add_edge g 8 10;
+  G.add_edge g 9 1;
+  G.add_edge g 9 7;
+  G.add_edge g 9 10;
+  G.add_edge g 10 1;
+  G.add_edge g 10 2;
+  G.add_edge g 11 1;
+  G.add_edge g 11 3
+
+let pp_comma p () = Format.(pp_print_char p ','; pp_print_space p ())
+let pp_set pf s =
+  Format.(fprintf pf "@[<hv 4>{%a}@]"
+            (pp_print_list ~pp_sep:pp_comma pp_print_int)
+            (Dominator.S.elements s))
+
+let () =
+  let idom = Dominator.compute_idom g 1 in
+  let domtree = Dominator.idom_to_dom_tree g idom in
+  let s1 = Dominator.S.of_list
+      (Dominator.dom_tree_to_nontrivial_dom 1 domtree)
+  in
+  let s2 = Dominator.dom_tree_to_snontrivial_dom 1 domtree in
+  Format.(printf "@[<v>%a@]@." (pp_print_list pp_set) [s1; s2])
+
diff --git a/tests/test_saps.expected b/tests/test_saps.expected
new file mode 100644
index 0000000..03771f5
--- /dev/null
+++ b/tests/test_saps.expected
@@ -0,0 +1,5 @@
+{3, 5}
+{1, 2, 3, 4, 5, 6}
+{B, C, E, I}
+{1, 2, 3, 4, 6, 7, 9, 12, 13, 14}
+{2, 4, 7, 12, 13}
diff --git a/tests/test_saps.ml b/tests/test_saps.ml
new file mode 100644
index 0000000..4afe887
--- /dev/null
+++ b/tests/test_saps.ml
@@ -0,0 +1,197 @@
+
+(* Test file for Connectivity.sstrong_articulation_points *)
+
+open Graph
+
+module Int = struct
+  type t = int
+  let compare = compare
+  let hash = Hashtbl.hash
+  let equal = (=)
+  let default = 0
+  end
+
+module String = struct
+  type t = string
+  let compare = compare
+  let hash = Hashtbl.hash
+  let equal = (=)
+  let default = "X"
+  end
+
+module GI = Imperative.Digraph.Concrete(Int)
+module CI = Components.Connectivity(Builder.I(GI))
+
+module GS = Persistent.Digraph.ConcreteBidirectional(String)
+module CS = Components.BiConnectivity(GS)
+
+(* From Fig. 1, Italiano, Laura, and Santaroni, "Finding strong bridges and
+   strong articulation points in linear time", TCS 447(2012). *)
+let g1 = GI.create ()
+let () =
+  GI.add_edge g1 1 2;
+  GI.add_edge g1 1 3;
+  GI.add_edge g1 2 3;
+  GI.add_edge g1 3 1;
+  GI.add_edge g1 3 4;
+  GI.add_edge g1 3 6;
+  GI.add_edge g1 3 7;
+  GI.add_edge g1 4 5;
+  GI.add_edge g1 5 1;
+  GI.add_edge g1 5 2;
+  GI.add_edge g1 5 4;
+  GI.add_edge g1 5 7;
+  GI.add_edge g1 6 5;
+  GI.add_edge g1 6 7;
+  GI.add_edge g1 7 6;
+  GI.add_edge g1 7 5
+
+(* From Fig. 2, Italiano, Laura, and Santaroni, "Finding strong bridges and
+   strong articulation points in linear time", TCS 447(2012). *)
+let g2 = GI.create ()
+let () =
+  GI.add_edge g2 1 2;
+  GI.add_edge g2 2 3;
+  GI.add_edge g2 3 4;
+  GI.add_edge g2 4 5;
+  GI.add_edge g2 5 6;
+  GI.add_edge g2 6 1
+
+(* From Fig. 7, Italiano, Laura, and Santaroni, "Finding strong bridges and
+   strong articulation points in linear time", TCS 447(2012). *)
+let g3 =
+  List.fold_left (fun g (src, dst) -> GS.add_edge g src dst) GS.empty
+    [("A", "B");
+     ("B", "A");
+     ("B", "E");
+     ("C", "A");
+     ("C", "B");
+     ("C", "D");
+     ("C", "E");
+     ("D", "C");
+     ("D", "E");
+     ("E", "C");
+     ("E", "D");
+     ("E", "J");
+     ("E", "I");
+     ("F", "D");
+     ("F", "J");
+     ("H", "F");
+     ("H", "I");
+     ("I", "H");
+     ("I", "J");
+     ("I", "E");
+     ("J", "E");
+     ("J", "F");
+     ("J", "I") ]
+
+(* From Fig. 3, Jaberi, "On computing the 2-vertex-connected components of
+   directed graphs", Discrete Applied Mathematicds 204(2016). *)
+let g4 = GI.create ()
+let () =
+  GI.add_edge g4 1 2;
+  GI.add_edge g4 2 3;
+  GI.add_edge g4 3 5;
+  GI.add_edge g4 3 9;
+  GI.add_edge g4 3 11;
+  GI.add_edge g4 3 12;
+  GI.add_edge g4 4 1;
+  GI.add_edge g4 4 5;
+  GI.add_edge g4 4 7;
+  GI.add_edge g4 4 9;
+  GI.add_edge g4 5 3;
+  GI.add_edge g4 5 4;
+  GI.add_edge g4 6 7;
+  GI.add_edge g4 6 9;
+  GI.add_edge g4 7 6;
+  GI.add_edge g4 7 8;
+  GI.add_edge g4 7 9;
+  GI.add_edge g4 8 6;
+  GI.add_edge g4 9 3;
+  GI.add_edge g4 9 4;
+  GI.add_edge g4 9 6;
+  GI.add_edge g4 9 7;
+  GI.add_edge g4 10 3;
+  GI.add_edge g4 10 11;
+  GI.add_edge g4 11 10;
+  GI.add_edge g4 11 12;
+  GI.add_edge g4 12 3;
+  GI.add_edge g4 12 10;
+  GI.add_edge g4 12 13;
+  GI.add_edge g4 13 14;
+  GI.add_edge g4 14 3;
+  GI.add_edge g4 14 10
+
+(* From Fig. 4, Jaberi, "On computing the 2-vertex-connected components of
+   directed graphs", Discrete Applied Mathematicds 204(2016). *)
+let g5 = GI.create ()
+let () =
+  GI.add_edge g5 0 1;
+  GI.add_edge g5 0 2;
+  GI.add_edge g5 1 0;
+  GI.add_edge g5 1 2;
+  GI.add_edge g5 2 0;
+  GI.add_edge g5 2 1;
+  GI.add_edge g5 2 3;
+  GI.add_edge g5 2 4;
+  GI.add_edge g5 2 5;
+  GI.add_edge g5 3 2;
+  GI.add_edge g5 3 4;
+  GI.add_edge g5 3 5;
+  GI.add_edge g5 4 2;
+  GI.add_edge g5 4 3;
+  GI.add_edge g5 4 5;
+  GI.add_edge g5 4 6;
+  GI.add_edge g5 4 7;
+  GI.add_edge g5 5 2;
+  GI.add_edge g5 5 3;
+  GI.add_edge g5 5 4;
+  GI.add_edge g5 6 4;
+  GI.add_edge g5 6 7;
+  GI.add_edge g5 7 4;
+  GI.add_edge g5 7 6;
+  GI.add_edge g5 7 8;
+  GI.add_edge g5 7 9;
+  GI.add_edge g5 7 11;
+  GI.add_edge g5 7 14;
+  GI.add_edge g5 8 7;
+  GI.add_edge g5 8 9;
+  GI.add_edge g5 8 10;
+  GI.add_edge g5 9 8;
+  GI.add_edge g5 9 10;
+  GI.add_edge g5 9 11;
+  GI.add_edge g5 9 13;
+  GI.add_edge g5 10 7;
+  GI.add_edge g5 10 11;
+  GI.add_edge g5 10 9;
+  GI.add_edge g5 11 7;
+  GI.add_edge g5 11 8;
+  GI.add_edge g5 11 10;
+  GI.add_edge g5 12 4;
+  GI.add_edge g5 12 13;
+  GI.add_edge g5 13 12;
+  GI.add_edge g5 13 14;
+  GI.add_edge g5 13 15;
+  GI.add_edge g5 14 2;
+  GI.add_edge g5 14 12;
+  GI.add_edge g5 15 12
+
+let pp_comma p () = Format.(pp_print_char p ','; pp_print_space p ())
+let pp_set pp_ele pf s =
+  Format.(fprintf pf "@[<hv 4>{%a}@]"
+            (pp_print_list ~pp_sep:pp_comma pp_ele) s)
+
+let () =
+  let saps1 = CI.sstrong_articulation_points g1 in
+  let saps2 = CI.sstrong_articulation_points g2 in
+  let saps3 = CS.sstrong_articulation_points g3 in
+  let saps4 = CI.sstrong_articulation_points g4 in
+  let saps5 = CI.sstrong_articulation_points g5 in
+  Format.(printf "@[<v>%a@,%a@,%a@,%a@,%a@]@."
+            (pp_set pp_print_int) (CI.S.elements saps1)
+            (pp_set pp_print_int) (CI.S.elements saps2)
+            (pp_set pp_print_string) (CS.S.elements saps3)
+            (pp_set pp_print_int) (CI.S.elements saps4)
+            (pp_set pp_print_int) (CI.S.elements saps5)
+         )
+
diff --git a/tests/test_topsort.ml b/tests/test_topsort.ml
index f74ad23..409a44e 100644
--- a/tests/test_topsort.ml
+++ b/tests/test_topsort.ml
@@ -66,7 +66,7 @@ let () = tests Topological.iter
 let rec pow a = function
   | 0 -> 1
   | 1 -> a
-  | n -> 
+  | n ->
     let b = pow a (n / 2) in
     b * b * (if n mod 2 = 0 then 1 else a)
 
@@ -80,4 +80,3 @@ let () =
     el := [n-1,0]; for i = 0 to n-2 do el := (i,i+1) :: !el done;
     test ~check:false Topological.iter n !el
   done
-