diff --git a/Project.toml b/Project.toml
index d6d4ae2..73bbb4f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -13,7 +13,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 [compat]
 BranchAndPrune = "0.2.1"
 ForwardDiff = "0.10, 1"
-IntervalArithmetic = "0.22.13, 0.23, 1"
+IntervalArithmetic = "1.0.3"
 Reexport = "1"
 StaticArrays = "1"
 julia = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 67fabeb..fb2b148 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,8 +1,10 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
-[sources.IntervalRootFinding]
-path = ".."
+[sources]
+IntervalRootFinding = {path = ".."}
diff --git a/docs/src/decorations.md b/docs/src/decorations.md
index bf718ee..1d7d83e 100644
--- a/docs/src/decorations.md
+++ b/docs/src/decorations.md
@@ -40,7 +40,7 @@ For example, `0.1` is famously not parsed as `0.1`,
 as `0.1` can not be represented exactly as a binary number
 (just like `1/3` can not be represented exactly as a decimal number).
 
-```julia-repl
+```jldoctest
 julia> big(0.1)
 0.1000000000000000055511151231257827021181583404541015625
 ```
diff --git a/docs/src/index.md b/docs/src/index.md
index 37fdd00..e2216ef 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,6 +1,6 @@
 ```@meta
 DocTestSetup = quote
-    using IntervalArithmetic, IntervalArithmetic.Symbols, IntervalRootFinding
+    using IntervalArithmetic, IntervalArithmetic.Symbols, IntervalRootFinding, StaticArrays
 end
 ```
 # `IntervalRootFinding.jl`
@@ -29,20 +29,40 @@ julia> using IntervalArithmetic, IntervalArithmetic.Symbols, IntervalRootFinding
 julia> rts = roots(x -> x^2 - 2x, 0..10)
 2-element Vector{Root{Interval{Float64}}}:
  Root([0.0, 3.73848e-8]_com_NG, :unknown)
+    Not converged: region size smaller than the tolerance
  Root([2.0, 2.0]_com_NG, :unique)
 ```
 
 The roots are returned as `Root` objects, containing an interval and the status of that interval, represented as a `Symbol`. There are two possible types of root status, as shown in the example:
   - `:unique`: the given interval contains *exactly one* root of the function,
   - `:unknown`: the given interval may or may not contain one or more roots; the algorithm used was unable to come to a conclusion.
+  In this case, additional information is stored in the root,
+  describing why the interval was not further bisected.
 
 The second status is still informative, since all regions of the original search interval *not* contained in *any* of the returned root intervals is guaranteed *not* to contain any root of the function. In the above example, we know that the function has no root in the interval $[2.1, 10]$, for example.
 
-There are several known situations where the uniqueness (and existence) of a solution cannot be determined by the interval algorithms used in the package:
+Most of the time, a root has status unknown because the convergence
+limits have been reached.
+This information is stored in the `convergence` field of the root,
+and can be either of
+  - `:max_iterartion`: the maximal number of iteration was reached
+  (`max_iteration` keyword);
+  - `:tolerance`: the interval is smaller than the specificed 
+  tolerance (`abstol` and `reltol` keywords,
+  for absolute and relative tolerance respectively).
+
+However, there are several known situations where the uniqueness
+(and existence) of a solution can never be determined
+by the interval algorithms used in the package:
+  - If the function used error at for specific intervals
+    (for example if comparison operators like `==` and `<` are used);
   - If the solution is on the boundary of the interval (as in the previous example);
   - If the derivative of the solution is zero at the solution.
 
-In particular, the second condition means that multiple roots cannot be proven to be unique. For example:
+The latest case happens for example when the function
+has multiple degenerate root,
+and then the unicity of the solution cannot be established.
+For example:
 
 ```jldoctest
 julia> g(x) = (x^2 - 2)^2 * (x^2 - 3)
@@ -52,7 +72,9 @@ julia> roots(g, -10..10)
 4-element Vector{Root{Interval{Float64}}}:
  Root([-1.73205, -1.73205]_com_NG, :unique)
  Root([-1.41421, -1.41421]_com, :unknown)
+    Not converged: region size smaller than the tolerance
  Root([1.41421, 1.41421]_com, :unknown)
+    Not converged: region size smaller than the tolerance
  Root([1.73205, 1.73205]_com_NG, :unique)
 ```
 
@@ -67,7 +89,7 @@ The initial search region is an array of interval.
 
 Here we give a 3D example:
 
-```julia-repl
+```jldoctest
 julia> function g( (x1, x2, x3) )
           return [
               x1^2 + x2^2 + x3^2 - 1,
@@ -82,28 +104,31 @@ julia> X = -5..5
 
 julia> rts = roots(g, [X, X, X])
 4-element Vector{Root{Vector{Interval{Float64}}}}:
- Root(Interval{Float64}[[-0.440763, -0.440762]_com_NG, [-0.866026, -0.866025]_com_NG, [0.236067, 0.236069]_com_NG], :unique)
- Root(Interval{Float64}[[-0.440763, -0.440762]_com_NG, [0.866025, 0.866026]_com_NG, [0.236067, 0.236069]_com_NG], :unique)
- Root(Interval{Float64}[[0.440762, 0.440763]_com_NG, [-0.866026, -0.866025]_com_NG, [0.236067, 0.236069]_com_NG], :unique)
- Root(Interval{Float64}[[0.440762, 0.440763]_com_NG, [0.866025, 0.866026]_com_NG, [0.236067, 0.236069]_com_NG], :unique)
+ Root(Interval{Float64}[[-0.440763, -0.440763]_com_NG, [-0.866025, -0.866025]_com_NG, [0.236068, 0.236068]_com_NG], :unique)
+ Root(Interval{Float64}[[-0.440763, -0.440763]_com_NG, [0.866025, 0.866025]_com_NG, [0.236068, 0.236068]_com_NG], :unique)
+ Root(Interval{Float64}[[0.440763, 0.440763]_com_NG, [-0.866025, -0.866025]_com_NG, [0.236068, 0.236068]_com_NG], :unique)
+ Root(Interval{Float64}[[0.440763, 0.440763]_com_NG, [0.866025, 0.866025]_com_NG, [0.236068, 0.236068]_com_NG], :unique)
 ```
 
 Thus, the system admits four unique roots in the box $[-5, 5]^3$.
 
 Moreover, the package is compatible with `StaticArrays.jl`.
 Usage of static arrays is recommended to increase performance.
-```julia-repl
+```jldoctest
 julia> using StaticArrays
 
 julia> h((x, y)) = SVector(x^2 - 4, y^2 - 16)
 h (generic function with 1 method)
 
+julia> X = -5..5
+[-5.0, 5.0]_com
+
 julia> roots(h, SVector(X, X))
 4-element Vector{Root{SVector{2, Interval{Float64}}}}:
- Root(Interval{Float64}[[-2.00001, -1.999999]_com_NG, [-4.00001, -3.999999]_com_NG], :unique)
- Root(Interval{Float64}[[-2.00001, -1.999999]_com_NG, [3.999999, 4.00001]_com_NG], :unique)
- Root(Interval{Float64}[[1.999999, 2.00001]_com_NG, [-4.00001, -3.999999]_com_NG], :unique)
- Root(Interval{Float64}[[1.999999, 2.00001]_com_NG, [3.999999, 4.00001]_com_NG], :unique)
+ Root(Interval{Float64}[[-2.0, -2.0]_com_NG, [-4.0, -4.0]_com_NG], :unique)
+ Root(Interval{Float64}[[-2.0, -2.0]_com_NG, [4.0, 4.0]_com_NG], :unique)
+ Root(Interval{Float64}[[2.0, 2.0]_com_NG, [-4.0, -4.0]_com_NG], :unique)
+ Root(Interval{Float64}[[2.0, 2.0]_com_NG, [4.0, 4.0]_com_NG], :unique)
 ```
 
 ### Vector types
@@ -167,14 +192,35 @@ julia> f( (x, y) ) = sin(x) * sin(y)
 f (generic function with 1 method)
 
 julia> ∇f(x) = gradient(f, x)  # gradient operator from the package
-(::#53) (generic function with 1 method)
+∇f (generic function with 1 method)
 
 julia> rts = roots(∇f, SVector(interval(-5, 6), interval(-5, 6)) ; abstol = 1e-5)
 25-element Vector{Root{SVector{2, Interval{Float64}}}}:
- Root(Interval{Float64}[[-4.7124, -4.71238]_com, [-4.7124, -4.71238]_com], :unique)
- Root(Interval{Float64}[[-3.1416, -3.14159]_com, [-3.1416, -3.14159]_com], :unique)
- ⋮
- [output skipped for brevity]
+ Root(Interval{Float64}[[-4.71239, -4.71239]_com, [-4.71239, -4.71239]_com], :unique)
+ Root(Interval{Float64}[[-3.14159, -3.14159]_com, [-3.14159, -3.14159]_com], :unique)
+ Root(Interval{Float64}[[-4.71239, -4.71239]_com, [-1.5708, -1.5708]_com], :unique)
+ Root(Interval{Float64}[[-3.14159, -3.14159]_com, [-3.88602e-10, 6.62844e-11]_com], :unique)
+ Root(Interval{Float64}[[-1.5708, -1.5708]_com, [-4.71239, -4.71239]_com], :unique)
+ Root(Interval{Float64}[[-3.74621e-10, 7.17246e-11]_com, [-3.14159, -3.14159]_com], :unique)
+ Root(Interval{Float64}[[-1.5708, -1.5708]_com, [-1.5708, -1.5708]_com], :unique)
+ Root(Interval{Float64}[[-6.04841e-12, 1.04386e-12]_com, [-6.04841e-12, 1.03449e-12]_com], :unique)
+ Root(Interval{Float64}[[-4.71239, -4.71239]_com, [1.5708, 1.5708]_com], :unique)
+ Root(Interval{Float64}[[-3.14159, -3.14159]_com, [3.14159, 3.14159]_com], :unique)
+ Root(Interval{Float64}[[-1.5708, -1.5708]_com, [1.5708, 1.5708]_com], :unique)
+ Root(Interval{Float64}[[-1.69679e-6, 3.8872e-7]_com, [3.14159, 3.14159]_com], :unique)
+ Root(Interval{Float64}[[-4.71239, -4.71239]_com, [4.71239, 4.71239]_com], :unique)
+ Root(Interval{Float64}[[-1.5708, -1.57079]_com, [4.71239, 4.71239]_com], :unique)
+ Root(Interval{Float64}[[1.5708, 1.5708]_com, [-4.71239, -4.71239]_com], :unique)
+ Root(Interval{Float64}[[3.14159, 3.14159]_com, [-3.14159, -3.14159]_com], :unique)
+ Root(Interval{Float64}[[1.5708, 1.5708]_com, [-1.5708, -1.5708]_com], :unique)
+ Root(Interval{Float64}[[3.14159, 3.14159]_com, [-3.5929e-14, 5.92849e-15]_com], :unique)
+ Root(Interval{Float64}[[4.71239, 4.71239]_com, [-4.71239, -4.71239]_com], :unique)
+ Root(Interval{Float64}[[4.71239, 4.71239]_com, [-1.5708, -1.57079]_com], :unique)
+ Root(Interval{Float64}[[1.5708, 1.5708]_com, [1.5708, 1.5708]_com], :unique)
+ Root(Interval{Float64}[[3.14159, 3.14159]_com, [3.14159, 3.14159]_com], :unique)
+ Root(Interval{Float64}[[1.57079, 1.5708]_com, [4.71239, 4.71239]_com], :unique)
+ Root(Interval{Float64}[[4.71239, 4.71239]_com, [1.57079, 1.5708]_com], :unique)
+ Root(Interval{Float64}[[4.71239, 4.71239]_com, [4.71239, 4.71239]_com], :unique)
 ```
 
 Now let's find the midpoints and plot them:
diff --git a/docs/src/internals.md b/docs/src/internals.md
index de8a29f..82bcad2 100644
--- a/docs/src/internals.md
+++ b/docs/src/internals.md
@@ -1,3 +1,8 @@
+```@meta
+DocTestSetup = quote
+    using IntervalArithmetic, IntervalArithmetic.Symbols, IntervalRootFinding
+end
+```
 # Internals
 
 This section describes some of the internal mechanism of the package and several ways to use them to customize a search.
@@ -37,14 +42,14 @@ The contractors are various methods to guarantee and refine the
 status of a root.
 The available contractors are `Bisection`, `Newton` or `Krawczyk`.
 
-```julia-repl
-julia> using IntervalArithmetic.Symbols
+```jldoctest
+julia> using IntervalRootFinding: contract
 
 julia> contract(Newton, sin, cos, Root(pi ± 0.001, :unknown))
-Root([3.14159, 3.1416]_com, :unique)
+Root([3.14159, 3.14159]_com, :unique)
 
 julia> contract(Newton, sin, cos, Root(2 ± 0.001, :unknown))
-Root([1.99899, 2.00101]_com, :empty)
+Root([1.999, 2.001]_com, :empty)
 ```
 
 ## RootProblem and search object
@@ -58,11 +63,20 @@ for example to log information at each iteration.
 For example, the following stops the search after 15 iterations and
 shows the state of the search at that point.
 
-```julia-repl
+```jldoctest
 julia> f(x) = exp(x) - sin(x)
 f (generic function with 1 method)
 
 julia> problem = RootProblem(f, interval(-10, 10))
+RootProblem
+  Contractor: Newton
+  Function: f
+  Search region: [-10.0, 10.0]_com
+  Search order: BranchAndPrune.BreadthFirst
+  Absolute tolerance: 1.0e-7
+  Relative tolerance: 0.0
+  Maximum iterations: 100000
+  Ignored errors: DataType[IntervalArithmetic.InconclusiveBooleanOperation]
 
 julia> state = nothing   # stores current state of the search
 
@@ -76,12 +90,16 @@ Branching
 └─ Branching
    ├─ Branching
    │  ├─ Branching
-   │  │  ├─ (:working, Root([-10.0, -8.7886]_com, :unknown))
-   │  │  └─ (:working, Root([-8.78861, -7.55813]_com, :unknown))
-   │  └─ (:final, Root([-6.28132, -6.28131]_com, :unique))
+   │  │  ├─ (:working, Root([-10.0, -8.78861]_com, :unknown)
+   │  │  │      Not converged: the root is still being processed)
+   │  │  └─ (:working, Root([-8.78861, -7.55814]_com, :unknown)
+   │  │         Not converged: the root is still being processed)
+   │  └─ (:final, Root([-6.28131, -6.28131]_com, :unique))
    └─ Branching
-      ├─ (:working, Root([-5.07783, -3.84734]_com, :unknown))
-      └─ (:working, Root([-3.84736, -2.5975]_com, :unknown))
+      ├─ (:working, Root([-5.07782, -3.84735]_com, :unknown)
+      │      Not converged: the root is still being processed)
+      └─ (:working, Root([-3.84735, -2.5975]_com, :unknown)
+             Not converged: the root is still being processed)
 ```
 
 The elements of the iteration are `SearchState` from the `BranchAndPrune.jl`
diff --git a/docs/src/roots.md b/docs/src/roots.md
index 9d24832..e55e57b 100644
--- a/docs/src/roots.md
+++ b/docs/src/roots.md
@@ -1,3 +1,8 @@
+```@meta
+DocTestSetup = quote
+    using IntervalArithmetic, IntervalArithmetic.Symbols, IntervalRootFinding
+end
+```
 # `roots` interface
 
 ## Methods
@@ -11,18 +16,19 @@ Both the Newton and Krawczyk methods can determine if a root is unique in an int
 
 The method used is given using the `contractor` keyword argument:
 
-```julia-repl
+```jldoctest
 julia> roots(log, -2..2 ; contractor = Newton)
 1-element Vector{Root{Interval{Float64}}}:
- Root([0.999999, 1.00001]_com, :unique)
+ Root([1.0, 1.0]_com, :unique)
 
 julia> roots(log, -2..2 ; contractor = Krawczyk)
 1-element Vector{Root{Interval{Float64}}}:
- Root([0.999999, 1.00001]_com, :unique)
+ Root([1.0, 1.0]_com, :unique)
 
 julia> roots(log, -2..2 ; contractor = Bisection)
 1-element Vector{Root{Interval{Float64}}}:
- Root([0.999999, 1.00001]_com, :unknown)
+ Root([1.0, 1.0]_com, :unknown)
+    Not converged: region size smaller than the tolerance
 ```
 
 Note that as shown in the example, the `log` function does not complain about being given an interval going outside of its domain. While this may be surprising, this is the expected behavior and no root will ever be found outside the domain of a function.
@@ -31,14 +37,14 @@ Note that as shown in the example, the `log` function does not complain about be
 
 Newton and Krawczyk methods require the function to be differentiable, but the derivative is usually computed automatically using forward-mode automatic differentiation, provided by the `ForwardDiff.jl` package. It is however possible to provide the derivative explicitly using the `derivative` keyword argument:
 
-```julia-repl
+```jldoctest
 julia> roots(log, -2..2 ; contractor = Newton, derivative = x -> 1/x)
 1-element Vector{Root{Interval{Float64}}}:
- Root([0.999999, 1.00001]_com_NG, :unique)
+ Root([1.0, 1.0]_com_NG, :unique)
 
 julia> roots(log, -2..2 ; contractor = Krawczyk, derivative = x -> 1/x)
 1-element Vector{Root{Interval{Float64}}}:
- Root([0.999999, 1.00001]_com_NG, :unique)
+ Root([1.0, 1.0]_com_NG, :unique)
 ```
 
 When providing the derivative explicitly, the computation is expected to be slightly faster, but the precision of the result is unlikely to be affected.
@@ -47,50 +53,50 @@ When providing the derivative explicitly, the computation is expected to be slig
 julia> using BenchmarkTools
 
 julia> @btime roots(log, -2..2 ; derivative = x -> 1/x)
-  4.814 μs (53 allocations: 3.16 KiB)
+  7.050 μs (129 allocations: 9.33 KiB)
 1-element Vector{Root{Interval{Float64}}}:
- Root([0.999999, 1.00001]_com_NG, :unique)
+ Root([1.0, 1.0]_com_NG, :unique)
 
 julia> @btime roots(log, -2..2)
-  5.767 μs (53 allocations: 3.16 KiB)
+  7.743 μs (129 allocations: 9.33 KiB)
 1-element Vector{Root{Interval{Float64}}}:
- Root([0.999999, 1.00001]_com, :unique)
+ Root([1.0, 1.0]_com, :unique)
 ```
 
 This may be useful in some special cases where `ForwardDiff.jl` is unable to compute the derivative of a function. Examples are complex functions and functions whose interval extension must be manually defined (e.g. special functions like `zeta`).
 
 In dimension greater than one, the derivative can be given as a function returning the Jacobi matrix:
 
-```julia-repl
+```jldoctest
 julia> f( (x, y) ) = [sin(x), cos(y)]
 f (generic function with 1 method)
 
 julia> df( (x, y) ) = [cos(x) 0 ; 0 -sin(y)]
+df (generic function with 1 method)
 
 julia> roots(f, [-3..3, -3..3] ; derivative = df)
 2-element Vector{Root{Vector{Interval{Float64}}}}:
- Root(Interval{Float64}[[-1.24409e-21, 1.0588e-22]_com_NG, [-1.5708, -1.57079]_com_NG], :unique)
- Root(Interval{Float64}[[-1.24409e-21, 1.0588e-22]_com_NG, [1.57079, 1.5708]_com_NG], :unique)
+ Root(Interval{Float64}[[-1.28378e-21, 1.58819e-22]_com_NG, [-1.5708, -1.5708]_com_NG], :unique)
+ Root(Interval{Float64}[[-1.28378e-21, 1.58819e-22]_com_NG, [1.5708, 1.5708]_com_NG], :unique)
 ```
 
 ## Tolerance
 
 An absolute tolerance for the search may be specified as the `abstol` keyword argument.
-Currently a method must first be provided in order to be able to choose the tolerance.
 
-```julia-repl
+```jldoctest
 julia> g(x) = sin(exp(x))
 g (generic function with 1 method)
 
 julia> roots(g, 0..2)
 2-element Vector{Root{Interval{Float64}}}:
- Root([1.14472, 1.14474]_com, :unique)
- Root([1.83787, 1.83788]_com, :unique)
+ Root([1.14473, 1.14473]_com, :unique)
+ Root([1.83788, 1.83788]_com, :unique)
 
 julia> roots(g, 0..2 ; abstol = 1e-1)
 2-element Vector{Root{Interval{Float64}}}:
  Root([1.14173, 1.15244]_com, :unique)
- Root([1.78757, 1.84273]_com, :unique)
+ Root([1.78757, 1.84272]_com, :unique)
 ```
 
 A lower tolerance may greatly reduce the computation time, at the cost of an increased number of returned roots having `:unknown` status:
@@ -100,31 +106,38 @@ julia> h(x) = cos(x) * sin(1 / x)
 h (generic function with 1 method)
 
 julia> @btime roots(h, 0.05..1)
-  79.500 μs (301 allocations: 13.97 KiB)
+  484.488 μs (12218 allocations: 590.76 KiB)
 6-element Vector{Root{Interval{Float64}}}:
  Root([0.0530516, 0.0530517]_com_NG, :unique)
- Root([0.0636619, 0.063662]_com_NG, :unique)
- Root([0.0795774, 0.0795775]_com_NG, :unique)
- Root([0.106103, 0.106104]_com_NG, :unique)
- Root([0.159154, 0.159156]_com_NG, :unique)
- Root([0.318309, 0.31831]_com_NG, :unique)
+ Root([0.063662, 0.063662]_com_NG, :unique)
+ Root([0.0795775, 0.0795775]_com_NG, :unique)
+ Root([0.106103, 0.106103]_com_NG, :unique)
+ Root([0.159155, 0.159155]_com_NG, :unique)
+ Root([0.31831, 0.31831]_com_NG, :unique)
 
 julia> @btime roots(h, 0.05..1 ; abstol = 1e-2)
-  48.500 μs (265 allocations: 11.61 KiB)
+  249.720 μs (6141 allocations: 297.80 KiB)
 6-element Vector{Root{Interval{Float64}}}:
- Root([0.0514445, 0.0531087]_com_NG, :unique)
- Root([0.0570253, 0.0641615]_com, :unknown)
+ Root([0.0514446, 0.0531086]_com_NG, :unique)
+ Root([0.0570254, 0.0641614]_com, :unknown)
+    Not converged: region size smaller than the tolerance
  Root([0.0785458, 0.0797797]_com_NG, :unknown)
- Root([0.104754, 0.107542]_com_NG, :unknown)
- Root([0.157236, 0.165989]_com_NG, :unknown)
- Root([0.31716, 0.319318]_com_NG, :unique)
+    Not converged: region size smaller than the tolerance
+ Root([0.104755, 0.107541]_com_NG, :unknown)
+    Not converged: region size smaller than the tolerance
+ Root([0.157236, 0.165988]_com_NG, :unknown)
+    Not converged: region size smaller than the tolerance
+ Root([0.317161, 0.319318]_com_NG, :unique)
 
 julia> @btime roots(h, 0.05..1 ; abstol = 1e-1)
-  17.300 μs (114 allocations: 5.16 KiB)
+  66.330 μs (1402 allocations: 69.10 KiB)
 3-element Vector{Root{Interval{Float64}}}:
- Root([0.04999999, 0.107542]_com, :unknown)
- Root([0.107541, 0.165989]_com, :unknown)
- Root([0.283803, 0.356099]_com_NG, :unknown)
+ Root([0.05, 0.107541]_com, :unknown)
+    Not converged: region size smaller than the tolerance
+ Root([0.107541, 0.165988]_com, :unknown)
+    Not converged: region size smaller than the tolerance
+ Root([0.283804, 0.356099]_com_NG, :unknown)
+    Not converged: region size smaller than the tolerance
 ```
 
 The last example shows a case where the tolerance was too large to be able to isolate the roots in distinct regions.
@@ -134,3 +147,40 @@ The last example shows a case where the tolerance was too large to be able to is
     For a root `x` of some function, if the absolute tolerance is smaller than `eps(x)` i.e. if `tol + x == x`, `roots` may never be able to converge to the required tolerance and the function may get stuck in an infinite loop.
     To avoid that, either increase `abstol` or `reltol`,
     or set a maximum number of iterations with the `max_iteration` keyword.
+
+## Error during computation
+
+By default, the `roots` function will ignore `IntervalArithmetic.InconclusiveBooleanOperation`
+errors and continue bisecting the given interval,
+hoping to eventually find regions that are small enough to avoid the error.
+
+This is useful when using code that is using comparisons,
+which are ill-defined for intervals
+(see [the IntervalArithmetic.jl documentation](https://juliaintervals.github.io/IntervalArithmetic.jl/stable/manual/usage/#Comparisons)
+for more information).
+
+For example, the following will find all roots,
+and flag the erroring singularity:
+```jldoctest bisect_on_error
+julia> f(x) = x > 2 ? 7 - x : x
+f (generic function with 1 method)
+
+julia> roots(f, -10 .. 10)
+3-element Vector{Root{Interval{Float64}}}:
+ Root([0.0, 0.0]_com, :unique)
+ Root([2.0, 2.0]_com, :unknown)
+    Not converged: region size smaller than the tolerance
+    Warning: an error was encountered during computation
+ Root([7.0, 7.0]_com_NG, :unique)
+```
+
+This behavior can be controlled by changing the `ignored_errors` field.
+
+To deactivate it, change it to an empty list,
+so that the process is interrupted as soon as an error is encountered.
+This is useful while debugging your code.
+
+```julia-repl
+julia> roots(f, -10 .. 10 ; ignored_errors = [])
+ERROR: ArgumentError: `<` is purposely not supported for overlapping intervals. See instead `strictprecedes`
+```
\ No newline at end of file
diff --git a/src/contractors.jl b/src/contractors.jl
index 56ac497..21032b5 100644
--- a/src/contractors.jl
+++ b/src/contractors.jl
@@ -3,6 +3,7 @@ function image_contains_zero(f, R::Root)
     R.status == :empty && return Root(X, :empty)
 
     imX = f(X)
+    size(imX) != size(X) && throw(DimensionMismatch("input and output dimensions of f must be the same"))
 
     if !(all(in_interval.(0, imX))) || isempty_region(imX)
         return Root(X, :empty)
@@ -57,7 +58,7 @@ function contract(::Type{Krawczyk}, f, derivative, X::AbstractVector)
     return mm - Y*f(mm) + (interval(I) - Y*J) * (X - mm)
 end
 
-function contract(::Type{C}, f, derivative, R::Root) where {C <: AbstractContractor}
+function contract(::Type{C}, f::F, derivative::Fp, R::Root) where {C <: AbstractContractor, F, Fp}
     # We first check with the simple bisection method
     # If we can prove it is empty at this point, we don't go further
     R2 = image_contains_zero(f, R)
@@ -105,5 +106,5 @@ function refine(root_problem::RootProblem{C}, R::Root) where C
         X = NX
     end
 
-    return Root(X, :unique)
+    return Root(X, :unique, :converged)
 end
\ No newline at end of file
diff --git a/src/root_object.jl b/src/root_object.jl
index a019157..8072c36 100644
--- a/src/root_object.jl
+++ b/src/root_object.jl
@@ -14,10 +14,23 @@ the `roots` function.
         representing an interval box) searched for roots.
   - `status`: the status of the region, valid values are `:empty`, `unknown`
         and `:unique`.
+  - `convergence`: the convergence status of the region. It is always `:converged`
+        for roots with status `:unique`,
+        and can be either `:max_iterartion` or `:tolerance` for roots with status `:unknown`,
+        depending on whether they stopped being processing due to reaching
+        the maximum number of iteration or the tolerance, respectively.
+  - `errored`: whether an error was encounter but ignored during the processing of this region.
+        The ignored errors are controlled by the `ignored_errors` field of the `RootProblem`.
 """
 struct Root{T}
     region::T
     status::Symbol
+    convergence::Symbol
+    errored::Bool
+end
+
+function Root(region, status::Symbol, convergence = :none, errored = false)
+    return Root(region, status, convergence, errored)
 end
 
 """
@@ -41,7 +54,24 @@ Return whether a `Root` is unique.
 """
 isunique(rt::Root{T}) where {T} = (rt.status == :unique)
 
-show(io::IO, rt::Root) = print(io, "Root($(rt.region), :$(rt.status))")
+function show(io::IO, rt::Root)
+    print(io, "Root($(rt.region), :$(rt.status))")
+    if rt.status == :unknown
+        if rt.convergence == :tolerance
+            print(io, "\n    Not converged: region size smaller than the tolerance")
+        elseif rt.convergence == :max_iterartion 
+            print(io, "\n    Not converged: reached maximal number of iterations")
+        elseif rt.convergence == :none
+            print(io, "\n    Not converged: the root is still being processed")
+        else
+            print(io, "\n    Not converged: unknown reason $(rt.convergence)")
+        end
+
+        if rt.errored
+            print(io, "\n    Warning: an error was encountered during computation")
+        end
+    end
+end
 
 ⊆(a::Interval, b::Root) = a ⊆ b.region
 ⊆(a::Root, b::Root) = a.region ⊆ b.region
diff --git a/src/roots.jl b/src/roots.jl
index baa44b2..a0e5b68 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -11,6 +11,7 @@ struct RootProblem{C, F, G, R, S, T}
     reltol::T
     max_iteration::Int
     where_bisect::T
+    ignored_errors::Vector{DataType}
 end
 
 """
@@ -44,11 +45,15 @@ Parameters
     of the remaining regions are smaller than `reltol * mag(interval)`.
     Default: `0.0`.
 - `max_iteration`: The maximum number of iteration, which also corresponds to
-    the maximum number of bisections allowed. Default: `typemax(Int)`.
+    the maximum number of bisections allowed. Default: `100_000`.
 - `where_bisect`: Value used to bisect the region. It is used to avoid
     bisecting exactly on zero when starting with symmetrical regions,
     often leading to having a solution directly on the boundary of a region,
     which prevent the contractor to prove it's unicity. Default: `127/256`.
+-`ignored_errors`: List of exceptions that are ignored during the processing
+    of a region. If the error is encoutered, it is discarded and the region is bisected
+    further.
+    Default: `[IntervalArithmetic.InconclusiveBooleanOperation]`.
 """
 RootProblem(f, region ; kwargs...) = RootProblem(f, Root(region, :unkown) ; kwargs...)
 
@@ -59,8 +64,9 @@ function RootProblem(
         search_order = BreadthFirst,
         abstol = 1e-7,
         reltol = 0.0,
-        max_iteration = typemax(Int),
-        where_bisect = 0.49609375)  # 127//256
+        max_iteration = 100_000,
+        where_bisect = 0.49609375,  # 127//256
+        ignored_errors = [IntervalArithmetic.InconclusiveBooleanOperation])
     
     N = length(root_region(root))
     if isnothing(derivative)
@@ -80,10 +86,23 @@ function RootProblem(
         abstol,
         reltol,
         max_iteration,
-        where_bisect
+        where_bisect,
+        convert(Vector{DataType}, ignored_errors)
     )
 end
 
+Base.show(io::IO, pb::RootProblem) = print(io, """
+    RootProblem
+      Contractor: $(pb.contractor)
+      Function: $(pb.f)
+      Search region: $(root_region(pb.region))
+      Search order: $(pb.search_order)
+      Absolute tolerance: $(pb.abstol)
+      Relative tolerance: $(pb.reltol)
+      Maximum iterations: $(pb.max_iteration)
+      Ignored errors: $(pb.ignored_errors)"""
+)
+
 function Base.iterate(root_problem::RootProblem, state = nothing)
     if isnothing(state)
         search = root_search(root_problem)
@@ -110,7 +129,14 @@ function under_tolerance(root_problem, root::Root)
 end
 
 function process(root_problem, root::Root)
-    contracted = contract(root_problem, root)
+    contracted = nothing
+    try
+        contracted = contract(root_problem, root)
+    catch err
+        !any(isa(err, Err) for Err in root_problem.ignored_errors) && rethrow()
+        contracted = Root(root.region, :unknown, :none, true)
+    end
+
     status = root_status(contracted)
 
     if status == :unique
@@ -122,8 +148,13 @@ function process(root_problem, root::Root)
 
     if status == :unknown
         # Avoid infinite division of intervals with singularity
-        istrivial(contracted.region) && under_tolerance(root_problem, root) && return :store, root
-        under_tolerance(root_problem, contracted) && return :store, contracted
+        if istrivial(contracted.region) && under_tolerance(root_problem, root)
+            return :store, Root(root.region, :unknown, :tolerance)
+        end
+
+        if under_tolerance(root_problem, contracted)
+            return :store, Root(contracted.region, :unknown, :tolerance, contracted.errored)
+        end
         
         return :branch, root
     else
@@ -173,7 +204,13 @@ function roots(f, region ; kwargs...)
         endstate.tree
     )
 
-    return vcat(result.final_regions, result.unfinished_regions)
+    rts = vcat(result.final_regions, result.unfinished_regions)
+    return map(rts) do rt
+        if rt.status == :unknown && rt.convergence == :none
+            return Root(rt.region, rt.status, :max_iterartion, rt.errored)
+        end
+        return rt
+    end
 end
 
 # Acting on complex `Interval`
diff --git a/test/roots.jl b/test/roots.jl
index 9246e47..9b3ba1c 100644
--- a/test/roots.jl
+++ b/test/roots.jl
@@ -296,7 +296,7 @@ end
 
     for abstol in abstols
         for reltol in reltols
-            rts = roots(f, interval(0, 1) ; abstol, reltol)
+            rts = roots(f, interval(0, 1) ; abstol, reltol, max_iteration = typemax(Int))
             regions = [rt.region for rt in rts if root_status(rt) == :unknown]
 
             @test all(
@@ -371,4 +371,24 @@ end
 
         @test all(isequal_interval.(myroots, unroots))
     end
+end
+
+@testset "Bisect on error" begin
+    f(x) = (x < 1 ? x : x^2)
+    
+    rts = roots(f, interval(-10, 10))
+    @test length(rts) == 2
+    @test rts[1].status == :unique
+    @test rts[1].convergence == :converged
+    @test !rts[1].errored
+
+    @test rts[2].status == :unknown
+    @test rts[2].convergence == :tolerance
+    @test rts[2].errored
+
+    @test_throws IntervalArithmetic.InconclusiveBooleanOperation roots(f, interval(-10, 10) ; ignored_errors = [])
+
+    g(x) = (x < 1 ? x : error())
+    rts = roots(g, interval(-10, 10) ; ignored_errors = [ErrorException, IntervalArithmetic.InconclusiveBooleanOperation])
+    @test any(isunique, rts)
 end
\ No newline at end of file