diff --git a/Project.toml b/Project.toml
index e5b52b60..16849fa2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,11 +27,11 @@ PolyhedraRecipesBaseExt = "RecipesBase"
 [compat]
 GenericLinearAlgebra = "0.2, 0.3"
 GeometryBasics = "0.2, 0.3, 0.4, 0.5"
-JuMP = "0.23, 1"
+JuMP = "1"
 LinearAlgebra = "1.6"
 MathOptInterface = "1"
 MutableArithmetics = "1"
-RecipesBase = "0.7, 0.8, 1.0"
+RecipesBase = "1.0"
 SparseArrays = "1.6"
-StaticArrays = "0.12, 1.0"
-julia = "1.6"
+StaticArrays = "1.0"
+julia = "1.10"
diff --git a/ext/PolyhedraGeometryBasicsExt.jl b/ext/PolyhedraGeometryBasicsExt.jl
index 7c4f0319..2a11a4bf 100644
--- a/ext/PolyhedraGeometryBasicsExt.jl
+++ b/ext/PolyhedraGeometryBasicsExt.jl
@@ -90,7 +90,7 @@ function _isdup(zray, triangles)
 end
 _isdup(poly, hidx, triangles) = _isdup(get(poly, hidx).a, triangles)
 
-function decompose_plane!(triangles::Vector, d::FullDim, zray, incident_points, incident_lines, incident_rays, exit_point::Function, counterclockwise::Function, rotate::Function)
+function decompose_plane!(triangles::Vector, d::FullDim, zray, incident_points, incident_lines, incident_rays, exit_point::Function, counterclockwise::Function, rotate::Function; tol)
     # xray should be the rightmost ray
     xray = nothing
     # xray should be the leftmost ray
@@ -98,7 +98,7 @@ function decompose_plane!(triangles::Vector, d::FullDim, zray, incident_points,
     isapproxzero(zray) && return
 
     # Checking rays
-    hull, lines, rays = _planar_hull(d, incident_points, incident_lines, incident_rays, counterclockwise, rotate)
+    hull, lines, rays = _planar_hull(d, incident_points, incident_lines, incident_rays, counterclockwise, rotate; tol)
     if isempty(lines)
         if length(hull) + length(rays) < 3
             return
@@ -190,16 +190,16 @@ function fulldecompose(triangles::Vector{_Tri{N,T}}) where {N,T}
     return pts, faces, ns
 end
 
-function fulldecompose(poly_geom::Mesh, ::Type{T}) where T
+function fulldecompose(poly_geom::Mesh, ::Type{T}; tol = Polyhedra._default_tol(T)) where T
     poly = poly_geom.polyhedron
     exit_point = scene(poly_geom, T)
     triangles = _Tri{2,T}[]
     z = StaticArrays.SVector(zero(T), zero(T), one(T))
-    decompose_plane!(triangles, FullDim(poly), z, collect(points(poly)), lines(poly), rays(poly), exit_point, counterclockwise, rotate)
+    decompose_plane!(triangles, FullDim(poly), z, collect(points(poly)), lines(poly), rays(poly), exit_point, counterclockwise, rotate; tol)
     return fulldecompose(triangles)
 end
 
-function fulldecompose(poly_geom::Mesh{3}, ::Type{T}) where T
+function fulldecompose(poly_geom::Mesh{3}, ::Type{T}; tol = Polyhedra._default_tol(T)) where T
     poly = poly_geom.polyhedron
     exit_point = scene(poly_geom, T)
     triangles = _Tri{3,T}[]
@@ -207,7 +207,7 @@ function fulldecompose(poly_geom::Mesh{3}, ::Type{T}) where T
         zray = get(poly, hidx).a
         counterclockwise(a, b) = dot(cross(a, b), zray)
         rotate(r) = cross(zray, r)
-        decompose_plane!(triangles, FullDim(poly), zray, incidentpoints(poly, hidx), incidentlines(poly, hidx), incidentrays(poly, hidx), exit_point, counterclockwise, rotate)
+        decompose_plane!(triangles, FullDim(poly), zray, incidentpoints(poly, hidx), incidentlines(poly, hidx), incidentrays(poly, hidx), exit_point, counterclockwise, rotate; tol)
     end
     for hidx in eachindex(hyperplanes(poly))
         decompose_plane(hidx)
diff --git a/src/comp.jl b/src/comp.jl
index 7f6d6f3b..e7a3a53b 100644
--- a/src/comp.jl
+++ b/src/comp.jl
@@ -5,6 +5,12 @@ _default_tol(::Type{<:Union{Integer, Rational}}) = MA.Zero()
 # of created everytime
 _default_tol(::Type{T}) where {T<:AbstractFloat} = Base.rtoldefault(T)
 
+_default_tol(::Type{T}, ::Type{U}) where {T,U} = _default_tol(promote_type(T, U))
+
+# We should **not** define `isless(::Float64, ::MutableArithmetics.Zero)`.
+# When we get such `MethodError`, it a nice error that always indicate some underlying issue
+# so we should fix that underlying issue, not add this method that would just hide bugs.
+
 isapproxzero(x::T; kws...) where {T<:Real} = x == zero(T)
 isapproxzero(x::T; tol=Base.rtoldefault(T)) where {T<:AbstractFloat} = abs(x) < tol
 isapproxzero(x::AbstractVector{T}; kws...) where {T<:Real} = isapproxzero(maximum(abs, x); kws...)
diff --git a/src/elements.jl b/src/elements.jl
index 1527a89e..becd9b16 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -280,22 +280,41 @@ for ElemT in [:HalfSpace, :HyperPlane, :Ray, :Line]
     end
 end
 
-ininterior(r::Ray, h::HalfSpace) = _neg(_dot(h.a, r))
-ininterior(l::Line, h::HalfSpace) = _neg(_dot(h.a, l))
-ininterior(p::AbstractVector, h::HalfSpace) = _lt(_dot(h.a, p), h.β)
+function ininterior(r::Ray{T}, h::HalfSpace{U}; tol = _default_tol(T, U)) where {T,U}
+    return _neg(_dot(h.a, r); tol)
+end
+
+function ininterior(l::Line{T}, h::HalfSpace{U}; tol = _default_tol(T, U)) where {T,U}
+    return _neg(_dot(h.a, l); tol)
+end
+
+function ininterior(p::AbstractVector{T}, h::HalfSpace{U}; tol = _default_tol(T, U)) where {T,U}
+    return _lt(_dot(h.a, p), h.β; tol)
+end
+
+function inrelativeinterior(p::VRepElement, h::HalfSpace; tol...)
+    return ininterior(p, h; tol...)
+end
+
+function Base.in(r::Ray{T}, h::HalfSpace{U}; tol = _default_tol(T, U)) where {T,U}
+    return _nonpos(_dot(h.a, r); tol)
+end
 
-inrelativeinterior(p::VRepElement, h::HalfSpace) = ininterior(p, h)
+function Base.in(l::Line{T}, h::HalfSpace{U}; tol = _default_tol(T, U)) where {T,U}
+    return _nonpos(_dot(h.a, l); tol)
+end
+
+function Base.in(p::AbstractVector{T}, h::HalfSpace{U}; tol = _default_tol(T, U)) where {T,U}
+    return _leq(_dot(h.a, p), h.β; tol)
+end
 
-Base.in(r::Ray, h::HalfSpace; kws...) = _nonpos(_dot(h.a, r); kws...)
-Base.in(l::Line, h::HalfSpace; kws...) = _nonpos(_dot(h.a, l); kws...)
-Base.in(p::AbstractVector, h::HalfSpace; kws...) = _leq(_dot(h.a, p), h.β; kws...)
+ininterior(p::VRepElement, h::HyperPlane; tol...) = false
 
-ininterior(p::VRepElement, h::HyperPlane) = false
-inrelativeinterior(p::VRepElement, h::HyperPlane) = p in h
+inrelativeinterior(p::VRepElement, h::HyperPlane; tol...) = return in(p, h; tol...)
 
-Base.in(r::Ray, h::HyperPlane; kws...) = isapproxzero(_dot(h.a, r); kws...)
-Base.in(l::Line, h::HyperPlane; kws...) = isapproxzero(_dot(h.a, l); kws...)
-Base.in(p::AbstractVector, h::HyperPlane; kws...) = _isapprox(_dot(h.a, p), h.β; kws...)
+Base.in(r::Ray, h::HyperPlane; tol...) = isapproxzero(_dot(h.a, r); tol...)
+Base.in(l::Line, h::HyperPlane; tol...) = isapproxzero(_dot(h.a, l); tol...)
+Base.in(p::AbstractVector, h::HyperPlane; tol...) = _isapprox(_dot(h.a, p), h.β; tol...)
 
 #function Base.vec(x::FixedVector{T}) where {T}
 #    y = Vector{T}(N)
diff --git a/src/iterators.jl b/src/iterators.jl
index bb6f0a93..674d1806 100644
--- a/src/iterators.jl
+++ b/src/iterators.jl
@@ -149,14 +149,14 @@ for (isVrep, elt, loop_singular) in [(true, :AbstractVector, :point),
 
         Returns the list of $plural incident to idx for the polyhedron `p`.
         """
-        $inc(p::Polyhedron{T}, idx) where {T} = get(p, IncidentElements{T, $elemtype(p)}(p, idx))
+        $inc(p::Polyhedron{T}, idx; tol = _default_tol(T)) where {T} = get(p, IncidentElements{T, $elemtype(p)}(p, idx); tol)
 
         """
             incident$(singular)indices(p::Polyhedron, idx)
 
-        Returns the list of the indices of $plural incident to idx for the polyhedron `p`.
+        Returns the list of the indices of $plural incident to `idx` for the polyhedron `p`.
         """
-        $incidx(p::Polyhedron{T}, idx; tol) where {T} = get(p, IncidentIndices{T, $elemtype(p)}(p, idx); tol)
+        $incidx(p::Polyhedron{T}, idx; tol = _default_tol(T)) where {T} = get(p, IncidentIndices{T, $elemtype(p)}(p, idx); tol)
 
         """
             $lenp($horvrep::$HorVRep)
diff --git a/src/polyhedron.jl b/src/polyhedron.jl
index d7c013a8..678b47f2 100644
--- a/src/polyhedron.jl
+++ b/src/polyhedron.jl
@@ -96,8 +96,8 @@ function setvrep! end
 Returns the `fulldim(p)`-dimensional hyper-volume of the polyhedron `p`.
 Returns `Inf` or `-one(T)` if it is infinite depending on whether the type `T` has an infinite value.
 """
-function volume(p::Polyhedron{T}) where T
-    return reduce(+, unscaled_volume_simplex(Δ) for Δ in triangulation(p);
+function volume(p::Polyhedron{T}; tol = _default_tol(T)) where T
+    return reduce(+, unscaled_volume_simplex(Δ) for Δ in triangulation(p; tol);
                   init=zero(T)) / factorial(fulldim(p))
 end
 
diff --git a/src/repelemop.jl b/src/repelemop.jl
index 7ad05466..6914ee71 100644
--- a/src/repelemop.jl
+++ b/src/repelemop.jl
@@ -57,21 +57,22 @@ end
 ######
 # IN #
 ######
-function _vinh(v::VRepElement, it::ElemIt{<:HRepElement}, infun)
+function _vinh(v::VRepElement, it::ElemIt{<:HRepElement}, infun; tol)
     for h in it
-        if !infun(v, h)
+        if !infun(v, h; tol)
             return false
         end
     end
     return true
 end
-function _vinh(v::VRepElement, hr::HRep, infun)
+
+function _vinh(v::VRepElement, hr::HRep, infun; kws...)
     # The line below fails on Julia v0.7. On the simplextest,
     # map(...) returns (false,) and then all((false,)) returns true
     # because (false,)[1] returns true...
     # all(map(it -> _vinh(v, it, infun), hreps(hr)))
     for it in hreps(hr)
-        if !_vinh(v, it, infun)
+        if !_vinh(v, it, infun; kws...)
             return false
         end
     end
@@ -89,7 +90,9 @@ If `h` is an halfspace, it returns whether ``\\langle a, x \\rangle \\le \\beta`
 
 Returns whether `p` is in `h`, e.g. in all the hyperplanes and halfspaces supporting `h`.
 """
-Base.in(v::VRepElement, hr::HRep) = _vinh(v, hr, Base.in)
+function Base.in(v::VRepElement{T}, hr::HRep{U}; tol = _default_tol(T, U)) where {T,U}
+    _vinh(v, hr, Base.in; tol)
+end
 
 """
     ininterior(p::VRepElement, h::HRepElement)
@@ -102,20 +105,24 @@ If `h` is an halfspace ``\\langle a, x \\rangle \\leq \\beta``, it returns wheth
 
 Returns whether `p` is in the interior of `h`, e.g. in the interior of all the hyperplanes and halfspaces supporting `h`.
 """
-ininterior(v::VRepElement, hr::HRep) = _vinh(v, hr, ininterior)
+function ininterior(v::VRepElement{T}, hr::HRep{U}; tol = _default_tol(T, U)) where {T,U}
+    return _vinh(v, hr, ininterior; tol)
+end
 
 """
-    inrelativeinterior(p::VRepElement, h::HRepElement)
+    inrelativeinterior(p::VRepElement, h::HRepElement; tol)
 
 Returns whether `p` is in the relative interior of `h`.
 If `h` is an hyperplane, it is equivalent to `p in h` since the relative interior of an hyperplane is itself.
 If `h` is an halfspace, it is equivalent to `ininterior(p, h)`.
 
-    inrelativeinterior(p::VRepElement, h::HRep)
+    inrelativeinterior(p::VRepElement, h::HRep; tol)
 
 Returns whether `p` is in the relative interior of `h`, e.g. in the relative interior of all the hyperplanes and halfspaces supporting `h`.
 """
-inrelativeinterior(v::VRepElement, hr::HRep) = _vinh(v, hr, inrelativeinterior)
+function inrelativeinterior(v::VRepElement{T}, hr::HRep{U}; tol = _default_tol(T, U)) where {T,U}
+    return _vinh(v, hr, inrelativeinterior; tol)
+end
 
 structural_nonzero_indices(a::SparseArrays.SparseVector) = SparseArrays.nonzeroinds(a)
 structural_nonzero_indices(a::AbstractVector) = eachindex(a)
@@ -157,13 +164,15 @@ function support_function(h::AbstractVector, rep::Rep, solver=Polyhedra.linear_o
     end
 end
 
-function _hinv(h::HRepElement, vr::ElemIt{<:VRepElement})
-    return all(v -> in(v, h), vr)
+function _hinv(h::HRepElement, vr::ElemIt{<:VRepElement}; kws...)
+    return all(v -> in(v, h; kws...), vr)
 end
-function _hinv(h::HRepElement, vr::VRep)
-    return all(v -> _hinv(h, v), vreps(vr))
+
+function _hinv(h::HRepElement, vr::VRep; kws...)
+    return all(v -> _hinv(h, v; kws...), vreps(vr))
 end
-function _hinh(h::HalfSpace, hr::HRep, solver)
+
+function _hinh(h::HalfSpace, hr::HRep, solver; tol)
     # ⟨a, x⟩ ≦ β -> if β < max ⟨a, x⟩ then h is outside
     model = support_function_model(h.a, hr, solver)
     MOI.optimize!(model)
@@ -173,30 +182,43 @@ function _hinh(h::HalfSpace, hr::HRep, solver)
     elseif term == MOI.INFEASIBLE
         return true
     elseif term == MOI.OPTIMAL
-        return _leq(MOI.get(model, MOI.ObjectiveValue()), h.β)
+        obj = MOI.get(model, MOI.ObjectiveValue())
+        # If `h` and `hr` have rational representations,
+        # the default tolerance is `MA.Zero`.
+        # However, if the solver only supports `Float64` then
+        # `obj` is `Float64` so we'll have to use some tolerance
+        if obj isa AbstractFloat && tol isa MA.Zero
+            tol = _default_tol(typeof(obj))
+        end
+        return _leq(obj, h.β; tol)
     else
         error("Cannot determine whether the polyhedron is contained in the",
               " halfspace or not because the linear program terminated with",
               " status $term.")
     end
 end
-function _hinh(h::HyperPlane, hr::HRep, solver)
-    _hinh(halfspace(h), hr, solver) && _hinh(halfspace(-h), hr, solver)
+function _hinh(h::HyperPlane, hr::HRep, solver; kws...)
+    _hinh(halfspace(h), hr, solver; kws...) && _hinh(halfspace(-h), hr, solver; kws...)
 end
 
-Base.issubset(vr::VRepresentation, h::HRepElement) = _hinv(h, vr)
-Base.issubset(hr::HRepresentation, h::HRepElement) = _hinh(h, hr)
+function Base.issubset(vr::VRepresentation{T}, h::HRepElement{U}; tol = _default_tol(T, U)) where {T,U}
+    return _hinv(h, vr; tol)
+end
+
+function Base.issubset(hr::HRepresentation{T}, h::HRepElement{U}; tol = _default_tol(T, U)) where {T,U}
+    return _hinh(h, hr; tol)
+end
 
 """
     issubset(p::Rep, h::HRepElement)
 
 Returns whether `p` is a subset of `h`, i.e. whether `h` supports the polyhedron `p`.
 """
-function Base.issubset(p::Polyhedron, h::HRepElement, solver=Polyhedra.linear_objective_solver(p))
+function Base.issubset(p::Polyhedron{T}, h::HRepElement{U}, solver=Polyhedra.linear_objective_solver(p); tol = _default_tol(T, U)) where {T,U}
     if vrepiscomputed(p)
-        _hinv(h, p)
+        _hinv(h, p; tol)
     else
-        _hinh(h, p, solver)
+        _hinh(h, p, solver; tol)
     end
 end
 
diff --git a/src/triangulation.jl b/src/triangulation.jl
index 65ed7a41..00bc63cb 100644
--- a/src/triangulation.jl
+++ b/src/triangulation.jl
@@ -35,7 +35,8 @@ function _triangulation(Δs, Δ, v_idx, h_idx, incident_idx, is_weak_adjacent, c
         end
     end
 end
-function triangulation_indices(p::Polyhedron; tol)
+
+function triangulation_indices(p::Polyhedron{T}; tol = _default_tol(T)) where {T}
     hasrays(p) && error("Triangulation only supported for polytope.")
     v_idx = eachindex(points(p))
     h_idx = eachindex(halfspaces(p))
@@ -46,6 +47,7 @@ function triangulation_indices(p::Polyhedron; tol)
     _triangulation(Δs, Δ, v_idx, h_idx, incident_idx, is_weak_adjacent, fulldim(p))
     return unique(Δs)
 end
-function triangulation(p::Polyhedron; tol)
+
+function triangulation(p::Polyhedron{T}; tol = _default_tol(T)) where {T}
     return map(Δ -> vrep(get.(p, Δ)), triangulation_indices(p; tol))
 end
diff --git a/src/vrep_optimizer.jl b/src/vrep_optimizer.jl
index 27726afc..eccef539 100644
--- a/src/vrep_optimizer.jl
+++ b/src/vrep_optimizer.jl
@@ -94,7 +94,7 @@ function MOI.optimize!(lpm::VRepOptimizer{T}) where T
         lpm.solution = first(points(prob))
     else
         better(a, b) = (lpm.objective_sense == MOI.MAX_SENSE ? a > b : a < b)
-        _better(a, b) = (lpm.objective_sense == MOI.MAX_SENSE ? _gt(a, b; tol) : _lt(a, b; tol))
+        _better(a, b) = (lpm.objective_sense == MOI.MAX_SENSE ? _gt(a, b; lpm.tol) : _lt(a, b; lpm.tol))
         bestobjval = zero(T)
         lpm.solution = nothing
         for r in allrays(prob)
diff --git a/test/lp_to_polyhedra.jl b/test/lp_to_polyhedra.jl
index f0774959..8c4b36e1 100644
--- a/test/lp_to_polyhedra.jl
+++ b/test/lp_to_polyhedra.jl
@@ -68,9 +68,9 @@ function MOI.get(mock::MockOptimizer{T}, attr::MOI.ConstraintPrimal,
 end
 function MOI.get(mock::MockOptimizer, ::MOI.ObjectiveValue)
     if mock.status == MOI.OPTIMAL
-        return _dot(mock.objective_func, mock.solution) + mock.objective_constant
+        return Polyhedra._dot(mock.objective_func, mock.solution) + mock.objective_constant
     elseif mock.status == MOI.DUAL_INFEASIBLE
-        return _dot(mock.objective_func, mock.solution)
+        return Polyhedra._dot(mock.objective_func, mock.solution)
     else
         error("No objective value available when termination status is $(mock.status).")
     end