diff --git a/perf/runbenchmarks.jl b/perf/runbenchmarks.jl
index 75867a79..fb892afe 100644
--- a/perf/runbenchmarks.jl
+++ b/perf/runbenchmarks.jl
@@ -39,24 +39,33 @@ using StaticArrays
 #   samples:          9809
 #   evals/sample:     1
 
-let
-    h = hrep([HalfSpace([ 1.,  1], 1),
-              HalfSpace([ 1., -1], 0),
-              HalfSpace([-1.,  0], 0)])
+using BenchmarkTools
+using Polyhedra
+
+function vector(T)
+    h = hrep([HalfSpace(T[ 1,  1], T(1)),
+              HalfSpace(T[ 1, -1], T(0)),
+              HalfSpace(T[-1,  0], T(0))])
     display(@benchmark(doubledescription($h)))
 end
 
-let
-    h = hrep([HalfSpace((@SVector [ 1,  1]), 1),
-              HalfSpace((@SVector [ 1, -1]), 0),
-              HalfSpace((@SVector [-1,  0]), 0)])
+using StaticArrays
+function static(T)
+    h = hrep([HalfSpace((@SVector T[ 1,  1]), T(1)),
+              HalfSpace((@SVector T[ 1, -1]), T(0)),
+              HalfSpace((@SVector T[-1,  0]), T(0))])
     display(@benchmark(doubledescription($h)))
+    @profview_allocs for _ in 1:10000
+        doubledescription(h)
+    end
 end
 
-let
-    h = hrep([ 1  1
-               1 -1
-              -1  0],
-             [1., 0, 0])
+function matrix(T)
+    h = hrep(T[
+        1  1
+        1 -1
+       -1  0],
+        T[1, 0, 0],
+    )
     display(@benchmark(doubledescription($h)))
 end
diff --git a/src/Polyhedra.jl b/src/Polyhedra.jl
index bec47008..5841f56a 100644
--- a/src/Polyhedra.jl
+++ b/src/Polyhedra.jl
@@ -9,8 +9,7 @@ using LinearAlgebra
 #      like RowEchelon.jl (slow) or LU (faster).
 import GenericLinearAlgebra
 
-import MutableArithmetics
-const MA = MutableArithmetics
+import MutableArithmetics as MA
 
 import MathOptInterface as MOI
 import MathOptInterface.Utilities as MOIU
diff --git a/src/center.jl b/src/center.jl
index 3ac6b476..31b2414a 100644
--- a/src/center.jl
+++ b/src/center.jl
@@ -19,7 +19,7 @@ function maximum_radius_with_center(h::HRep{T}, center) where T
                     return zero(T)
                 end
             else
-                new_radius = (hs.β - center ⋅ hs.a) / n
+                new_radius = (hs.β - _dot(center, hs.a)) / n
                 if radius === nothing
                     radius = new_radius
                 else
diff --git a/src/comp.jl b/src/comp.jl
index 643459e3..7f6d6f3b 100644
--- a/src/comp.jl
+++ b/src/comp.jl
@@ -1,14 +1,21 @@
+# We use `Zero` so that it is allocation-free, while `zero(Rational{BigInt})` would be costly
+_default_tol(::Type{<:Union{Integer, Rational}}) = MA.Zero()
+# `rtoldefault(BigFloat)` is quite slow and allocates a lot so we hope that
+# `tol` will be called at some upper level function and passed in here instead
+# of created everytime
+_default_tol(::Type{T}) where {T<:AbstractFloat} = Base.rtoldefault(T)
+
 isapproxzero(x::T; kws...) where {T<:Real} = x == zero(T)
-isapproxzero(x::T; ztol=Base.rtoldefault(T)) where {T<:AbstractFloat} = abs(x) < ztol
-isapproxzero(x::AbstractVector{T}; kws...) where {T<:Real} = isapproxzero(maximum(abs.(x)); kws...)
+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...)
 #isapproxzero(x::Union{SymPoint, Ray, Line}; kws...) = isapproxzero(coord(x); kws...)
 isapproxzero(x::Union{Ray, Line}; kws...) = isapproxzero(coord(x); kws...)
 isapproxzero(h::HRepElement; kws...) = isapproxzero(h.a; kws...) && isapproxzero(h.β; kws...)
 
-_isapprox(x::Union{T, AbstractArray{T}}, y::Union{T, AbstractArray{T}}) where {T<:Union{Integer, Rational}} = x == y
+_isapprox(x::Union{T, AbstractArray{T}}, y::Union{T, AbstractArray{T}}; tol) where {T<:Union{Integer, Rational}} = x == y
 # I check with zero because isapprox(0, 1e-100) is false...
 # but isapprox(1e-100, 2e-100) should be false
-_isapprox(x, y) = (isapproxzero(x) ? isapproxzero(y) : (isapproxzero(y) ? isapproxzero(x) : isapprox(x, y)))
+_isapprox(x, y; tol) = (isapproxzero(x; tol) ? isapproxzero(y; tol) : (isapproxzero(y; tol) ? isapproxzero(x; tol) : isapprox(x, y; rtol = tol)))
 
 #Base.isapprox(r::SymPoint, s::SymPoint) = _isapprox(coord(r), coord(s)) || _isapprox(coord(r), -coord(s))
 function _scaleray(r::Union{Line, Ray}, s::Union{Line, Ray})
@@ -16,12 +23,12 @@ function _scaleray(r::Union{Line, Ray}, s::Union{Line, Ray})
     cs = coord(s)
     cr * sum(abs.(cs)), cs * sum(abs.(cr))
 end
-function Base.isapprox(r::Line, s::Line)
+function _isapprox(r::Line, s::Line; tol)
     rs, ss = _scaleray(r, s)
-    _isapprox(rs, ss) || _isapprox(rs, -ss)
+    _isapprox(rs, ss; tol) || _isapprox(rs, -ss; tol)
 end
-function Base.isapprox(r::Ray, s::Ray)
-    _isapprox(_scaleray(r, s)...)
+function _isapprox(r::Ray, s::Ray; tol)
+    _isapprox(_scaleray(r, s)...; tol)
 end
 # TODO check that Vec in GeometryTypes also does that
 function _scalehp(h1, h2)
@@ -33,28 +40,28 @@ function Base.:(==)(h1::HyperPlane, h2::HyperPlane)
     (a1, β1), (a2, β2) = _scalehp(h1, h2)
     (a1 == a2 && β1 == β2) || (a1 == -a2 && β1 == -β2)
 end
-function Base.isapprox(h1::HyperPlane, h2::HyperPlane)
+function _isapprox(h1::HyperPlane, h2::HyperPlane; tol)
     (a1, β1), (a2, β2) = _scalehp(h1, h2)
-    (_isapprox(a1, a2) && _isapprox(β1, β2)) || (_isapprox(a1, -a2) && _isapprox(β1, -β2))
+    (_isapprox(a1, a2; tol) && _isapprox(β1, β2; tol)) || (_isapprox(a1, -a2; tol) && _isapprox(β1, -β2; tol))
 end
 function Base.:(==)(h1::HalfSpace, h2::HalfSpace)
     (a1, β1), (a2, β2) = _scalehp(h1, h2)
     a1 == a2 && β1 == β2
 end
-function Base.isapprox(h1::HalfSpace, h2::HalfSpace)
+function _isapprox(h1::HalfSpace, h2::HalfSpace; tol)
     (a1, β1), (a2, β2) = _scalehp(h1, h2)
-    _isapprox(a1, a2) && _isapprox(β1, β2)
+    _isapprox(a1, a2; tol) && _isapprox(β1, β2; tol)
 end
 
-_lt(x::T, y::T) where {T<:Real} = x < y
-_lt(x::T, y::T) where {T<:AbstractFloat} = x < y && !_isapprox(x, y)
-_lt(x::S, y::T) where {S<:Real,T<:Real} = _lt(promote(x, y)...)
-_gt(x::S, y::T) where {S<:Real, T<:Real} = _lt(y, x)
-_leq(x::T, y::T) where {T<:Real} = x <= y
-_leq(x::T, y::T) where {T<:AbstractFloat} = x <= y || _isapprox(x, y)
-_leq(x::S, y::T) where {S<:Real,T<:Real} = _leq(promote(x, y)...)
-_geq(x::T, y::T) where {T<:Real} = _leq(y, x)
-_pos(x::T) where {T<:Real} = _gt(x, zero(T))
-_neg(x::T) where {T<:Real} = _lt(x, zero(T))
-_nonneg(x::T) where {T<:Real} = _geq(x, zero(T))
-_nonpos(x::T) where {T<:Real} = _leq(x, zero(T))
+_lt(x::T, y::T; tol) where {T<:Real} = x < y
+_lt(x::T, y::T; tol) where {T<:AbstractFloat} = x < y && !_isapprox(x, y; tol)
+_lt(x::S, y::T; tol) where {S<:Real,T<:Real} = _lt(promote(x, y)...; tol)
+_gt(x::S, y::T; tol) where {S<:Real, T<:Real} = _lt(y, x; tol)
+_leq(x::T, y::T; tol) where {T<:Real} = x <= y
+_leq(x::T, y::T; tol) where {T<:AbstractFloat} = x <= y || _isapprox(x, y; tol)
+_leq(x::S, y::T; tol) where {S<:Real,T<:Real} = _leq(promote(x, y)...; tol)
+_geq(x::T, y::T; tol) where {T<:Real} = _leq(y, x; tol)
+_pos(x::T; tol) where {T<:Real} = _gt(x, zero(T); tol)
+_neg(x::T; tol) where {T<:Real} = _lt(x, zero(T); tol)
+_nonneg(x::T; tol) where {T<:Real} = _geq(x, zero(T); tol)
+_nonpos(x::T; tol) where {T<:Real} = _leq(x, zero(T); tol)
diff --git a/src/doubledescription.jl b/src/doubledescription.jl
index 3eb7dbbb..e8ba5ea6 100644
--- a/src/doubledescription.jl
+++ b/src/doubledescription.jl
@@ -19,7 +19,7 @@ function dualfullspace(rep::Representation{T}) where T
 end
 
 """
-    doubledescription(h::HRepresentation)
+    doubledescription(h::HRepresentation; tol)
 
 Computes the V-representation of the polyhedron represented by `h` using the Double-Description algorithm [MRTT53, FP96].
 It maintains a list of the points, rays and lines that are in the resulting representation
@@ -305,21 +305,21 @@ function set_in!(data, I, index)
         end
     end
 end
-function add_element!(data, k, el, tight)
+function add_element!(data, k, el, tight; tol)
     cutoff = nothing
     for i in reverse(1:k)
         if data.cutline[i] !== nothing
             el = line_project(el, data.cutline[i], data.halfspaces[i])
             # The line is in all halfspaces from `i+1` up so projecting with it does not change it.
             push!(tight, i)
-            return add_adjacent_element!(data, i, k, el, data.lineray[i], tight)
+            return add_adjacent_element!(data, i, k, el, data.lineray[i], tight; tol)
         end
         # Could avoid computing `dot` twice between `el` and the halfspace here.
-        if !(el in data.halfspaces[i])
+        if !in(el, data.halfspaces[i]; tol)
             cutoff = i
             break
         end
-        if el in hyperplane(data.halfspaces[i])
+        if in(el, hyperplane(data.halfspaces[i]); tol)
             push!(tight, i)
         end
     end
@@ -327,22 +327,22 @@ function add_element!(data, k, el, tight)
     set_in!(data, (index.cutoff + 1):k, index)
     return index
 end
-function project_onto_affspace(data, offset, el, hyperplanes)
+function project_onto_affspace(data, offset, el, hyperplanes; tol)
     for i in reverse(eachindex(hyperplanes))
         line = data.cutline[offset + i]
         h = hyperplanes[i]
         if line !== nothing
             el = line_project(el, line, h)
-        elseif !(el in h)
+        elseif !in(el, h; tol)
             # e.g. 0x1 + 0x2 = 1 or -1 = x1 = 1, ...
             return nothing
         end
     end
     return el
 end
-function add_adjacent_element!(data, i, k, el, parent, tight)
-    index = add_element!(data, i - 1, el, tight)
-    add_intersection!(data, index, parent, nothing, i - 1)
+function add_adjacent_element!(data, i, k, el, parent, tight; tol)
+    index = add_element!(data, i - 1, el, tight; tol)
+    add_intersection!(data, index, parent, nothing, i - 1; tol)
     set_in!(data, i:k, index)
     return index
 end
@@ -377,11 +377,11 @@ Combine `el1` and `el2` which are on two different sides of the halfspace `h` so
 that the combination is in the corresponding hyperplane.
 """
 function combine(h::HalfSpace, el1::VRepElement, el2::VRepElement)
-    combine(h.β, el1, h.a ⋅ el1, el2, h.a ⋅ el2)
+    combine(h.β, el1, _dot(h.a, el1), el2, _dot(h.a, el2))
 end
 
 """
-    add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)
+    add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1; tol)
 
 Add the intersection of the elements of indices `idx1` and `idx2` that belong to
 the hyperplane corresponding the halfspace of index `hp_idx` (see 3.2 (ii) of
@@ -396,9 +396,9 @@ to avoid adding redundant elements.
 Double description method revisited
 *Combinatorics and computer science*, *Springer*, **1996**, 91-111
 """
-function add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)
+function add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1; tol)
     if idx1.cutoff > idx2.cutoff
-        return add_intersection!(data, idx2, idx1, hp_idx, hs_idx)
+        return add_intersection!(data, idx2, idx1, hp_idx, hs_idx; tol)
     end
     i = idx2.cutoff
     # Condition (c_k) in [FP96]
@@ -418,9 +418,13 @@ function add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)
     #      in case `rj, rj'` have inherited adjacency
     tight = tight_halfspace_indices(data, idx1) ∩ tight_halfspace_indices(data, idx2)
     push!(tight, i)
-    return add_adjacent_element!(data, i, i, newel, idx1, tight)
+    return add_adjacent_element!(data, i, i, newel, idx1, tight; tol)
 end
 
+# `MA.operate` is going to redirect to an implementation that
+# exploits the mutable API of `BigFloat` of `BigInt`
+_dot(a, b) = MA.operate(LinearAlgebra.dot, a, b)
+
 _shift(el::AbstractVector, line::Line) = el + Polyhedra.coord(line)
 _shift(el::Line, line::Line) = el + line
 _shift(el::Ray, line::Line) = el + Polyhedra.Ray(Polyhedra.coord(line))
@@ -428,21 +432,21 @@ function _λ_proj(r::VStruct, line::Line, h::HRepElement)
     # Line or ray `r`:
     # (r + λ * line) ⋅ h.a == 0
     # λ = -r ⋅ h.a / (line ⋅ h.a)
-    return -r ⋅ h.a / (line ⋅ h.a)
+    return -_dot(r, h.a) / _dot(line, h.a)
 end
 function _λ_proj(x::AbstractVector, line::Line, h::HRepElement)
     # Point `x`:
     # (x + λ * line) ⋅ h.a == h.β
     # λ = (h.β - x ⋅ h.a) / (line ⋅ h.a)
-    return (h.β - x ⋅ h.a) / (line ⋅ h.a)
+    return (h.β - _dot(x, h.a)) / _dot(line, h.a)
 end
 function line_project(el, line, h)
     λ = _λ_proj(el, line, h)
     return Polyhedra.simplify(_shift(el, λ * line))
 end
-function hline(data, line::Line, i, h)
-    value = h.a ⋅ line
-    if !Polyhedra.isapproxzero(value)
+function hline(data, line::Line, i, h; tol)
+    value = _dot(h.a, line)
+    if !Polyhedra.isapproxzero(value; tol)
         if data.cutline[i] === nothing
             # Make `lineray` point inward
             data.cutline[i] = value > 0 ? -line : line
@@ -456,7 +460,7 @@ function hline(data, line::Line, i, h)
 end
 
 # TODO remove solver arg `_`, it is kept to avoid breaking code
-function doubledescription(hr::HRepresentation, _ = nothing)
+function doubledescription(hr::HRepresentation{T}, solver = nothing; tol = _default_tol(T)) where {T}
     v = Polyhedra.dualfullspace(hr)
     hps = Polyhedra.lazy_collect(hyperplanes(hr))
     hss = Polyhedra.lazy_collect(halfspaces(hr))
@@ -464,12 +468,12 @@ function doubledescription(hr::HRepresentation, _ = nothing)
     for line in lines(v)
         cut = false
         for i in reverse(eachindex(hps))
-            cut, line = hline(data, line, nhalfspaces(hr) + i, hps[i])
+            cut, line = hline(data, line, nhalfspaces(hr) + i, hps[i]; tol)
             cut && break
         end
         if !cut
             for i in reverse(eachindex(hss))
-                cut, line = hline(data, line, i, hss[i])
+                cut, line = hline(data, line, i, hss[i]; tol)
                 cut && break
             end
         end
@@ -484,14 +488,14 @@ function doubledescription(hr::HRepresentation, _ = nothing)
         line = data.cutline[i]
         if line !== nothing
             ray = Polyhedra.Ray(Polyhedra.coord(line))
-            data.lineray[i] = add_element!(data, i - 1, ray, BitSet((i + 1):length(hss)))
+            data.lineray[i] = add_element!(data, i - 1, ray, BitSet((i + 1):length(hss)); tol)
         end
     end
     @assert isone(npoints(v))
     # Add the origin
-    orig = project_onto_affspace(data, nhalfspaces(hr), first(points(v)), hps)
+    orig = project_onto_affspace(data, nhalfspaces(hr), first(points(v)), hps; tol)
     if orig !== nothing
-        add_element!(data, nhalfspaces(hr), orig, resized_bitset(data))
+        add_element!(data, nhalfspaces(hr), orig, resized_bitset(data); tol)
     end
     for i in reverse(eachindex(hss))
         if data.cutline[i] === nothing && isempty(data.cutpoints[i]) && isempty(data.cutrays[i])
@@ -507,11 +511,11 @@ function doubledescription(hr::HRepresentation, _ = nothing)
             # Catches new adjacent rays, see 3.2 (ii) of [FP96]
             for p1 in data.pin[i], p2 in data.pin[i]
                 if p1.cutoff < p2.cutoff
-                    add_intersection!(data, p1, p2, i)
+                    add_intersection!(data, p1, p2, i; tol)
                 end
             end
             for p in data.pin[i], r in data.rin[i]
-                add_intersection!(data, p, r, i)
+                add_intersection!(data, p, r, i; tol)
             end
         end
         deleteat!(data.cutpoints, i)
@@ -521,7 +525,7 @@ function doubledescription(hr::HRepresentation, _ = nothing)
                 # We encounter both `r1, r2` and `r2, r1`.
                 # Break this symmetry with:
                 if r1.cutoff < r2.cutoff
-                    add_intersection!(data, r1, r2, i)
+                    add_intersection!(data, r1, r2, i; tol)
                 end
             end
         end
diff --git a/src/elements.jl b/src/elements.jl
index 31f2bef6..1527a89e 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -161,13 +161,17 @@ Base.:+(r::Line, s::Line) = Line(r.a + s.a)
 Base.:+(r::Ray, s::Ray) = Ray(r.a + s.a)
 Base.:+(p::AbstractVector, r::Ray) = p + coord(r)
 
-for op in [:dot, :cross]
-    @eval begin
-        LinearAlgebra.$op(x::VStruct, y) = $op(x.a, y)
-        LinearAlgebra.$op(x, y::VStruct) = $op(x, y.a)
-        LinearAlgebra.$op(x::VStruct, y::VStruct) = $op(x.a, y.a)
-    end
-end
+LinearAlgebra.cross(x::VStruct, y) = cross(x.a, y)
+LinearAlgebra.cross(x, y::VStruct) = cross(x, y.a)
+LinearAlgebra.cross(x::VStruct, y::VStruct) = cross(x.a, y.a)
+
+LinearAlgebra.dot(x::VStruct, y) = _dot(x, y)
+LinearAlgebra.dot(x, y::VStruct) = _dot(x, y)
+LinearAlgebra.dot(x::VStruct, y::VStruct) = _dot(x, y)
+
+_dot(x::VStruct, y) = _dot(x.a, y)
+_dot(x, y::VStruct) = _dot(x, y.a)
+_dot(x::VStruct, y::VStruct) = _dot(x.a, y.a)
 
 for ElT in (:HyperPlane, :HalfSpace, :Line, :Ray)
     @eval begin
@@ -276,22 +280,22 @@ for ElemT in [:HalfSpace, :HyperPlane, :Ray, :Line]
     end
 end
 
-ininterior(r::Ray, h::HalfSpace) = _neg(h.a ⋅ r)
-ininterior(l::Line, h::HalfSpace) = _neg(h.a ⋅ l)
-ininterior(p::AbstractVector, h::HalfSpace) = _lt(h.a ⋅ p, h.β)
+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.β)
 
 inrelativeinterior(p::VRepElement, h::HalfSpace) = ininterior(p, h)
 
-Base.in(r::Ray, h::HalfSpace) = _nonpos(h.a ⋅ r)
-Base.in(l::Line, h::HalfSpace) = _nonpos(h.a ⋅ l)
-Base.in(p::AbstractVector, h::HalfSpace) = _leq(h.a ⋅ p, h.β)
+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) = false
 inrelativeinterior(p::VRepElement, h::HyperPlane) = p in h
 
-Base.in(r::Ray, h::HyperPlane) = isapproxzero(h.a ⋅ r)
-Base.in(l::Line, h::HyperPlane) = isapproxzero(h.a ⋅ l)
-Base.in(p::AbstractVector, h::HyperPlane) = _isapprox(h.a ⋅ p, h.β)
+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...)
 
 #function Base.vec(x::FixedVector{T}) where {T}
 #    y = Vector{T}(N)
@@ -344,7 +348,7 @@ lift(v::AbstractVector{T}) where {T} = pushbefore(v, one(T))
 translate(p::AbstractVector, v) = p + v
 translate(r::VStruct, v) = r
 
-translate(h::ElemT, p) where {ElemT<:HRepElement} = ElemT(h.a, h.β + h.a ⋅ p)
+translate(h::ElemT, p) where {ElemT<:HRepElement} = ElemT(h.a, h.β + _dot(h.a, p))
 
 _simplify(a::AbstractVector) = a
 _simplify(a::AbstractVector, β) = a, β
diff --git a/src/incidence.jl b/src/incidence.jl
index 7c931f53..39bf3f5c 100644
--- a/src/incidence.jl
+++ b/src/incidence.jl
@@ -1,19 +1,19 @@
-isincident(v::VRepElement, h::HRepElement) = v in hyperplane(h)
-isincident(h::HRepElement, v::VRepElement) = v in hyperplane(h)
-isincident(p::Polyhedron, v::Index, h::Index) = isincident(get(p, v), get(p, h))
+isincident(v::VRepElement, h::HRepElement; tol) = in(v, hyperplane(h); tol)
+isincident(h::HRepElement, v::VRepElement; tol) = in(v, hyperplane(h); tol)
+isincident(p::Polyhedron, v::Index, h::Index; tol) = isincident(get(p, v), get(p, h); tol)
 
 abstract type Incident{T, ElemT<:RepElement{T}, PT<:Polyhedron{T}, IdxT<:Index{T}} end
 
-function Base.get(p::Polyhedron{T}, inc::Incident{T, ElemT}) where {T, ElemT}
+function Base.get(p::Polyhedron{T}, inc::Incident{T, ElemT}; tol) where {T, ElemT}
     el = get(p, inc.idx)
     incs = _inctype(inc)[]
     for idx in Indices{T, ElemT}(p)
         eli = get(p, idx)
-        if isincident(el, eli)
+        if isincident(el, eli; tol)
             push!(incs, _incel(inc, idx, eli))
         end
     end
-    incs
+    return incs
 end
 
 struct IncidentElements{T, ElemT, PT, IdxT} <: Incident{T, ElemT, PT, IdxT}
diff --git a/src/interval.jl b/src/interval.jl
index f31dbc23..b7e49152 100644
--- a/src/interval.jl
+++ b/src/interval.jl
@@ -17,8 +17,8 @@ mutable struct Interval{T, AT, D} <: Polyhedron{T}
     length::T
 end
 
-Interval{T, AT, D}(d::FullDim, it::HIt...) where {T, AT, D} = _hinterval(Intersection{T, AT, D}(d, it...), AT, D)
-Interval{T, AT, D}(d::FullDim, it::VIt...) where {T, AT, D} = _vinterval(Hull{T, AT, D}(d, it...), AT, D)
+Interval{T, AT, D}(d::FullDim, it::HIt...; tol = _default_tol(T)) where {T, AT, D} = _hinterval(Intersection{T, AT, D}(d, it...), AT, D; tol)
+Interval{T, AT, D}(d::FullDim, it::VIt...; tol = _default_tol(T)) where {T, AT, D} = _vinterval(Hull{T, AT, D}(d, it...), AT, D; tol)
 
 # If AT is an SVector, it will be StaticArrays.Size{(1,)}
 # otherwise, it will be 1
@@ -35,8 +35,8 @@ surface(::Interval{T}) where {T} = zero(T)
 volume(p::Interval) = p.length
 Base.isempty(p::Interval) = isempty(p.vrep)
 
-function _interval(AT, haslb::Bool, lb::T, hasub::Bool, ub::T, isempty::Bool) where {T}
-    if haslb && hasub && _gt(lb, ub)
+function _interval(AT, haslb::Bool, lb::T, hasub::Bool, ub::T, isempty::Bool; tol) where {T}
+    if haslb && hasub && _gt(lb, ub; tol)
         isempty = true
     end
     hps = HyperPlane{T, AT}[]
@@ -47,7 +47,7 @@ function _interval(AT, haslb::Bool, lb::T, hasub::Bool, ub::T, isempty::Bool) wh
     if !isempty
         if hasub
             push!(ps, StaticArrays.SVector(ub))
-            if haslb && _isapprox(lb, ub)
+            if haslb && _isapprox(lb, ub; tol)
                 push!(hps, HyperPlane(StaticArrays.SVector(one(T)), ub))
             else
                 push!(hss, HalfSpace(StaticArrays.SVector(one(T)), ub))
@@ -64,7 +64,7 @@ function _interval(AT, haslb::Bool, lb::T, hasub::Bool, ub::T, isempty::Bool) wh
             end
         end
         if haslb
-            if !_isapprox(lb, ub)
+            if !_isapprox(lb, ub; tol)
                 push!(hss, HalfSpace(StaticArrays.SVector(-one(T)), -lb))
                 push!(ps, StaticArrays.SVector(lb))
             end
@@ -84,7 +84,7 @@ function Interval{T, AT, D}(haslb::Bool, lb::T, hasub::Bool, ub::T, isempty::Boo
     return Interval{T, AT, D}(_interval(AT, haslb, lb, hasub, ub, isempty)...)
 end
 
-function _hinterval(rep::HRep{T}, ::Type{AT}) where {T, AT}
+function _hinterval(rep::HRep{T}, ::Type{AT}; tol) where {T, AT}
     haslb = false
     lb = zero(T)
     hasub = false
@@ -108,8 +108,8 @@ function _hinterval(rep::HRep{T}, ::Type{AT}) where {T, AT}
     end
     for hp in hyperplanes(rep)
         α = hp.a[1]
-        if isapproxzero(α)
-            if !isapproxzero(hp.β)
+        if isapproxzero(α; tol)
+            if !isapproxzero(hp.β; tol)
                 empty = true
             end
         else
@@ -119,7 +119,7 @@ function _hinterval(rep::HRep{T}, ::Type{AT}) where {T, AT}
     end
     for hs in halfspaces(rep)
         α = hs.a[1]
-        if isapproxzero(α)
+        if isapproxzero(α; tol)
             if hs.β < 0
                 empty = true
             end
@@ -129,14 +129,14 @@ function _hinterval(rep::HRep{T}, ::Type{AT}) where {T, AT}
             _setub(hs.β / α)
         end
     end
-    return _interval(AT, haslb, lb, hasub, ub, empty)
+    return _interval(AT, haslb, lb, hasub, ub, empty; tol)
 end
 
-function _hinterval(rep::HRep{T}, ::Type{AT}, D) where {T, AT}
-    return Interval{T, AT, D}(_hinterval(rep, AT)...)
+function _hinterval(rep::HRep{T}, ::Type{AT}, D; tol) where {T, AT}
+    return Interval{T, AT, D}(_hinterval(rep, AT; tol)...)
 end
 
-function _vinterval(v::VRep{T}, ::Type{AT}) where {T, AT}
+function _vinterval(v::VRep{T}, ::Type{AT}; tol) where {T, AT}
     haslb = true
     lb = zero(T)
     hasub = true
@@ -171,20 +171,20 @@ function _vinterval(v::VRep{T}, ::Type{AT}) where {T, AT}
             end
         end
     end
-    return _interval(AT, haslb, lb, hasub, ub, isempty)
+    return _interval(AT, haslb, lb, hasub, ub, isempty; tol)
 end
 
-function _vinterval(rep::VRep{T}, ::Type{AT}, D) where {T, AT}
-    return Interval{T, AT, D}(_vinterval(rep, AT)...)
+function _vinterval(rep::VRep{T}, ::Type{AT}, D; tol) where {T, AT}
+    return Interval{T, AT, D}(_vinterval(rep, AT; tol)...)
 end
 
-Interval{T, AT, D}(p::HRepresentation{T}) where {T, AT, D} = _hinterval(p, AT, D)
-Interval{T, AT, D}(p::VRepresentation{T}) where {T, AT, D} = _vinterval(p, AT, D)
-function Interval{T, AT, D}(p::Polyhedron{T}) where {T, AT, D}
+Interval{T, AT, D}(p::HRepresentation{T}; tol = _default_tol(T)) where {T, AT, D} = _hinterval(p, AT, D; tol)
+Interval{T, AT, D}(p::VRepresentation{T}; tol = _default_tol(T)) where {T, AT, D} = _vinterval(p, AT, D; tol)
+function Interval{T, AT, D}(p::Polyhedron{T}; tol = _default_tol(T)) where {T, AT, D}
     if hrepiscomputed(p)
-        _hinterval(p, AT, D)
+        _hinterval(p, AT, D; tol)
     else
-        _vinterval(p, AT, D)
+        _vinterval(p, AT, D; tol)
     end
 end
 
@@ -205,16 +205,16 @@ function removehredundancy!(::Interval, args...; kws...) end
 function removevredundancy!(::Interval, args...; kws...) end
 
 sethrep!(p::Interval, h::HRep) = resethrep!(p, h)
-function resethrep!(p::Interval{T, AT}, h::HRep) where {T, AT}
-    hnew, v, volume = _hinterval(h, AT)
+function resethrep!(p::Interval{T, AT}, h::HRep{U}; tol = _default_tol(U)) where {T, U, AT}
+    hnew, v, volume = _hinterval(h, AT; tol)
     p.hrep = hnew
     p.vrep = v
     p.length = volume
     return p
 end
 setvrep!(p::Interval, v::VRep) = resetvrep!(p, v)
-function resetvrep!(p::Interval{T, AT}, v::VRep) where {T, AT}
-    h, vnew, volume = _vinterval(v, AT)
+function resetvrep!(p::Interval{T, AT}, v::VRep{U}; tol = _default_tol(U)) where {T, U, AT}
+    h, vnew, volume = _vinterval(v, AT; tol)
     p.hrep = h
     p.vrep = vnew
     p.length = volume
diff --git a/src/iterators.jl b/src/iterators.jl
index 97627b30..bb6f0a93 100644
--- a/src/iterators.jl
+++ b/src/iterators.jl
@@ -156,7 +156,7 @@ for (isVrep, elt, loop_singular) in [(true, :AbstractVector, :point),
 
         Returns the list of the indices of $plural incident to idx for the polyhedron `p`.
         """
-        $incidx(p::Polyhedron{T}, idx) where {T} = get(p, IncidentIndices{T, $elemtype(p)}(p, idx))
+        $incidx(p::Polyhedron{T}, idx; tol) where {T} = get(p, IncidentIndices{T, $elemtype(p)}(p, idx); tol)
 
         """
             $lenp($horvrep::$HorVRep)
diff --git a/src/linearity.jl b/src/linearity.jl
index 3fcab565..8cbed30c 100644
--- a/src/linearity.jl
+++ b/src/linearity.jl
@@ -70,11 +70,11 @@ linearize(r::Ray) = line(r)
 _aff_push!(aff::HAffineSpace, h) = intersect!(aff, h)
 _aff_push!(aff::VLinearSpace, l) = convexhull!(aff, l)
 
-function _detect_opposite_element(aff, non_opposite, element; kws...)
+function _detect_opposite_element(aff, non_opposite, element; tol)
     element = remproj(element, aff)
-    if !isapproxzero(element; kws...) && !any(el -> el ≈ element, non_opposite)
+    if !isapproxzero(element; tol) && !any(el -> _isapprox(el, element; tol), non_opposite)
         lin = linearize(element)
-        i = findfirst(el -> linearize(el) ≈ lin, non_opposite)
+        i = findfirst(el -> _isapprox(linearize(el), lin; tol), non_opposite)
         if i === nothing
             push!(non_opposite, element)
         else
@@ -86,9 +86,9 @@ function _detect_opposite_element(aff, non_opposite, element; kws...)
         return false
     end
 end
-function _detect_opposite_elements(aff, non_opposite, elements; kws...)
+function _detect_opposite_elements(aff, non_opposite, elements; tol)
     newlin = true
-    for i in 1:fulldim(aff) # could use `while newlin` but `for`-loop is safer.
+    for _ in 1:fulldim(aff) # could use `while newlin` but `for`-loop is safer.
         newlin || break
         newlin = false
         empty!(non_opposite)
@@ -96,25 +96,25 @@ function _detect_opposite_elements(aff, non_opposite, elements; kws...)
         # Remove rays/halfspaces that become zero and detect new lines/hyperplanes
         # with rays/halfspaces that becomes opposite to each other.
         for element in elements
-            newlin |= _detect_opposite_element(aff, non_opposite, element; kws...)
+            newlin |= _detect_opposite_element(aff, non_opposite, element; tol)
         end
     end
 end
 
 """
-    detect_new_linearities(rep::HRepresentation{T}, solver; verbose=0, ztol=Base.rtoldefault(T)) where {T}
+    detect_new_linearities(rep::HRepresentation{T}, solver; verbose=0, tol=Base.rtoldefault(T)) where {T}
 
 Given a polyhedron with H-representation `rep`, detect whether a new hyperplane can be generated from the halfspaces in `halfspaces` using an linear program solved by `solver`.
 The method is similar to the method used for lines described as follows.
 This function is automatically called by `removehredundancy` if a solver is provided.
 
-    detect_new_linearities(rep::VRepresentation{T}, solver; verbose=0, ztol=Base.rtoldefault(T)) where {T}
+    detect_new_linearities(rep::VRepresentation{T}, solver; verbose=0, tol=Base.rtoldefault(T)) where {T}
 
 Given a cone defined by the V-representation `rep` (ignoring the points in the representation if any), detect whether a new line can be generated from the rays in `rays` using an linear program solved by `solver`.
 The method is as follows (suppose `lines` is empty for simplicity).
 This function is automatically called by `removevredundancy` if a solver is provided.
 
-The keyword argument `ztol` is used as a tolerance to decide whether a number is zero.
+The keyword argument `tol` is used as a tolerance to decide whether a number is zero.
 
 If there was a line `l` in the cone, it would mean that there exist `μ >= 0` and `ν >= 0` such that
 `Σ μ_i r_i = l` and `Σ ν_i r_i = -l`. We deduce from this that `Σ λ_i r_i = 0` where `λ = μ + ν`.
@@ -231,7 +231,7 @@ _linearity_space(h::HRepresentation, current) = affinehull(h, current)
 _linearity_space(v::VRepresentation, current) = linespace(v, current)
 
 struct OppositeMockOptimizer end
-function _detect_linearity(rep::Representation, solver; kws...)
+function _detect_linearity(rep::Representation, solver; verbose = 0, tol)
     aff = _linearity_space(rep, true)
     if _hasnonlinearity(rep)
         if solver === nothing
@@ -244,9 +244,9 @@ Set a solver if you believe that the polyhedron may have more linearity.
         end
         if solver == OppositeMockOptimizer
             els = _nonlin_type(rep)[]
-            _detect_opposite_elements(aff, els, _nonlinearity(rep); kws...)
+            _detect_opposite_elements(aff, els, _nonlinearity(rep); tol)
         else
-            new_lins = detect_new_linearities(rep, solver; kws...)
+            new_lins = detect_new_linearities(rep, solver; verbose, tol)
             if new_lins === nothing
                 empty!(aff)
                 els = eltype(_nonlinearity(rep))[]
@@ -273,8 +273,8 @@ Return a new H-representation with linearity detected using `solver`.
 
 The remaining keyword arguments `kws` are passed to [`detect_new_linearities`](@ref).
 """
-function detecthlinearity(hr::HRepresentation, solver; kws...)
-    aff, hs = _detect_linearity(hr, solver; kws...)
+function detecthlinearity(hr::HRepresentation{T}, solver; verbose = 0, tol = _default_tol(float(T))) where {T}
+    aff, hs = _detect_linearity(hr, solver; verbose = 0, tol)
     typeof(hr)(FullDim(hr), aff.hyperplanes, hs)
 end
 
@@ -285,7 +285,7 @@ Return a new V-representation with linearity detected using `solver`.
 
 The remaining keyword arguments `kws` are passed to [`detect_new_linearities`](@ref).
 """
-function detectvlinearity(vr::VRepresentation, solver; kws...)
-    aff, rays = _detect_linearity(vr, solver; kws...)
+function detectvlinearity(vr::VRepresentation{T}, solver; tol = _default_tol(float(T)), kws...) where {T}
+    aff, rays = _detect_linearity(vr, solver; tol, kws...)
     typeof(vr)(FullDim(vr), preps(vr)..., aff.lines, rays)
 end
diff --git a/src/planar.jl b/src/planar.jl
index 1a6e96ec..9ae45a86 100644
--- a/src/planar.jl
+++ b/src/planar.jl
@@ -92,7 +92,7 @@ end
 
 _colinear(counterclockwise, a, b, c) = isapproxzero(counterclockwise(b - a, c - a))
 
-function _planar_hull(d::FullDim, points, lines, rays, counterclockwise, rotate)
+function _planar_hull(d::FullDim, points, lines, rays, counterclockwise, rotate; tol)
     line = nothing
     lineleft = false
     lineright = false
@@ -123,10 +123,10 @@ function _planar_hull(d::FullDim, points, lines, rays, counterclockwise, rotate)
             if xray === nothing
                 xray = yray = coord(r)
             else
-                if Line(xray) ≈ linearize(r) && !(Ray(xray) ≈ r)
+                if _isapprox(Line(xray), linearize(r); tol) && !_isapprox(Ray(xray), r; tol)
                     line = Line(xray)
                     checkleftright(Ray(yray))
-                elseif Line(yray) ≈ linearize(r) && !(Ray(yray) ≈ r)
+                elseif _isapprox(Line(yray), linearize(r); tol) && !_isapprox(Ray(yray), r; tol)
                     line = Line(yray)
                     checkleftright(Ray(xray))
                 else
@@ -188,7 +188,7 @@ function _planar_hull(d::FullDim, points, lines, rays, counterclockwise, rotate)
         else
             append!(_points, half_points)
             push!(_rays, Ray(xray))
-            if !(Ray(xray) ≈ Ray(yray))
+            if !(_isapprox(Ray(xray), Ray(yray); tol))
                 push!(_rays, Ray(yray))
             end
         end
@@ -215,7 +215,7 @@ end
 counterclockwise(x, y) = x[1] * y[2] - x[2] * y[1]
 rotate(x) = convert(typeof(x), StaticArrays.SVector(-x[2], x[1]))
 
-function planar_hull(vr::VRepresentation)
+function planar_hull(vr::VRepresentation{T}; tol = _default_tol(T)) where {T}
     d = FullDim(vr)
-    vrep(_planar_hull(FullDim(vr), collect(points(vr)), lines(vr), rays(vr), counterclockwise, rotate)...; d = d)
+    vrep(_planar_hull(FullDim(vr), collect(points(vr)), lines(vr), rays(vr), counterclockwise, rotate; tol)...; d = d)
 end
diff --git a/src/polyhedron.jl b/src/polyhedron.jl
index ec2d180a..d7c013a8 100644
--- a/src/polyhedron.jl
+++ b/src/polyhedron.jl
@@ -119,12 +119,12 @@ end
 Returns the center of mass of `p`, represented as a `Vector{T}` of length `fulldim(p)`.
 Throws an error if `p` is degenerate.
 """
-function center_of_mass(p::Polyhedron{T}) where T
+function center_of_mass(p::Polyhedron{T}; tol = _default_tol(T)) where T
     # Implementation strategy: For a simplex, the center of mass coincides with
     # the centroid which is easy to compute.  So we triangulate `p` into
     # simplices and compute a volume-weighted average of these simplices'
     # centers of mass.
-    simplices = triangulation(p)
+    simplices = triangulation(p; tol)
     isempty(simplices) && error("Tried to compute center of mass of a degenerate polyhedron")
     unscaled_center = zeros(T, fulldim(p))
     unscaled_vol = zero(T)
diff --git a/src/redundancy.jl b/src/redundancy.jl
index 48eaa3cd..6855d93c 100644
--- a/src/redundancy.jl
+++ b/src/redundancy.jl
@@ -402,9 +402,9 @@ end
 # V-redundancy
 # If p is an H-representation, nl needs to be given otherwise if p is a Polyhedron, it can be asked to p.
 # TODO nlines should be the number of non-redundant lines so something similar to dim
-function isredundant(p::HRep{T}, v::Union{AbstractVector, Line, Ray}; strongly = false, nl::Int=nlines(p), solver=nothing) where {T}
+function isredundant(p::HRep{T}, v::Union{AbstractVector, Line, Ray}; strongly = false, nl::Int=nlines(p), solver=nothing, tol = _default_tol(T)) where {T}
     # v is in every hyperplane otherwise it would not be valid
-    hcount = nhyperplanes(p) + count(h -> v in hyperplane(h), halfspaces(p))
+    hcount = nhyperplanes(p) + count(h -> in(v, hyperplane(h); tol), halfspaces(p))
     strong = (isray(v) ? fulldim(p) - 1 : fulldim(p)) - nl
     bound = strongly ? min(strong, 1) : strong
     if hcount < bound
@@ -412,7 +412,7 @@ function isredundant(p::HRep{T}, v::Union{AbstractVector, Line, Ray}; strongly =
     else
         A = Matrix{T}(undef, hcount, fulldim(p))
         offset = coord_matrix!(A, hyperplanes(p), _ -> true, 0)
-        offset = coord_matrix!(A, halfspaces(p), h -> v in hyperplane(h), offset)
+        offset = coord_matrix!(A, halfspaces(p), h -> in(v, hyperplane(h); tol), offset)
         @assert offset == size(A, 1)
         return rank(A) < bound
     end
@@ -422,10 +422,10 @@ isredundant(p::HRep{T}, v::Line; strongly = false, nl::Int=nlines(p), solver=not
 
 # H-redundancy
 # If p is a V-representation, nl needs to be given otherwise if p is a Polyhedron, it can be asked to p.
-function isredundant(p::VRep{T}, h::HRepElement; strongly = false, d::Int=dim(p), solver=nothing) where {T}
+function isredundant(p::VRep{T}, h::HRepElement; strongly = false, d::Int=dim(p), solver=nothing, tol = _default_tol(T)) where {T}
     checkvconsistency(p)
     hp = hyperplane(h)
-    pcount = count(p -> p in hp, points(p))
+    pcount = count(p -> in(p, hp; tol), points(p))
     # every line is in h, otherwise it would not be valid
     rcount = nlines(p) + count(r -> r in hp, rays(p))
     if pcount < min(d, 1) || (!strongly && pcount + rcount < d)
@@ -437,7 +437,7 @@ function isredundant(p::VRep{T}, h::HRepElement; strongly = false, d::Int=dim(p)
         if pcount > 1
             orig = nothing
             for x in points(p)
-                if x in hp
+                if in(x, hp; tol)
                     if orig === nothing
                         orig = x
                     else
@@ -495,8 +495,8 @@ function removeduplicates(vrep::VPolytope)
     typeof(vrep)(FullDim(vrep), premovedups(vrep, emptyspace(vrep))...)
 end
 removeduplicates(vrep::VLinearSpace, solver = nothing) = detectvlinearity(vrep, solver)
-function removeduplicates(vrep::VRepresentation, solver = nothing)
-    aff, rs = _detect_linearity(vrep, solver)
+function removeduplicates(vrep::VRepresentation{T}, solver = nothing; tol = _default_tol(T)) where {T}
+    aff, rs = _detect_linearity(vrep, solver; tol)
     typeof(vrep)(FullDim(vrep), premovedups(vrep, aff)..., aff.lines, rs)
 end
 
diff --git a/src/repelemop.jl b/src/repelemop.jl
index 8049d900..7ad05466 100644
--- a/src/repelemop.jl
+++ b/src/repelemop.jl
@@ -282,9 +282,9 @@ function Base.intersect(v::VRepresentation{T}, h::HRepElement) where {T}
     #   \  In CDD and ConvexHull.jl, they only consider segments that
     #    \ belongs to the boundary which is way faster than what we do here
     for p in pins
-        ap = h.a ⋅ p
+        ap = _dot(h.a, p)
         for q in pout
-            aq = h.a ⋅ q
+            aq = _dot(h.a, q)
             λ = (aq - h.β) / (aq - ap)
             @assert 0 <= λ <= 1
             push!(pinp, simplify(λ * p + (1 - λ) * q))
@@ -292,9 +292,9 @@ function Base.intersect(v::VRepresentation{T}, h::HRepElement) where {T}
     end
     # Similar but with rays
     for r in rins
-        ar = h.a ⋅ r
+        ar = _dot(h.a, r)
         for s in rout
-            as = h.a ⋅ s
+            as = _dot(h.a, s)
             # should take
             # λ = as / (as - ar)
             @assert 0 <= as / (as - ar) <= 1
@@ -313,9 +313,9 @@ function Base.intersect(v::VRepresentation{T}, h::HRepElement) where {T}
     # Similar but with one point and one ray
     for (ps, rs) in ((pins, rout), (pout, rins))
         for p in ps
-            ap = h.a ⋅ p
+            ap = _dot(h.a, p)
             for r in rs
-                ar = h.a ⋅ r
+                ar = _dot(h.a, r)
                 λ = (h.β - ap) / ar
                 push!(pinp, p + λ * r)
             end
diff --git a/src/triangulation.jl b/src/triangulation.jl
index 1dddc5c5..65ed7a41 100644
--- a/src/triangulation.jl
+++ b/src/triangulation.jl
@@ -35,17 +35,17 @@ function _triangulation(Δs, Δ, v_idx, h_idx, incident_idx, is_weak_adjacent, c
         end
     end
 end
-function triangulation_indices(p::Polyhedron)
+function triangulation_indices(p::Polyhedron; tol)
     hasrays(p) && error("Triangulation only supported for polytope.")
     v_idx = eachindex(points(p))
     h_idx = eachindex(halfspaces(p))
     Δ = eltype(v_idx)[]
     Δs = typeof(Δ)[]
-    incident_idx = Dict(h => Set(incidentpointindices(p, h)) for h in h_idx)
+    incident_idx = Dict(h => Set(incidentpointindices(p, h; tol)) for h in h_idx)
     is_weak_adjacent = Dict((hi, hj) => !isempty(incident_idx[hi] ∩ incident_idx[hj]) for hi in h_idx for hj in h_idx)
     _triangulation(Δs, Δ, v_idx, h_idx, incident_idx, is_weak_adjacent, fulldim(p))
     return unique(Δs)
 end
-function triangulation(p::Polyhedron)
-    return map(Δ -> vrep(get.(p, Δ)), triangulation_indices(p))
+function triangulation(p::Polyhedron; tol)
+    return map(Δ -> vrep(get.(p, Δ)), triangulation_indices(p; tol))
 end
diff --git a/src/vrep_optimizer.jl b/src/vrep_optimizer.jl
index 730165ad..27726afc 100644
--- a/src/vrep_optimizer.jl
+++ b/src/vrep_optimizer.jl
@@ -23,10 +23,14 @@ mutable struct VRepOptimizer{T} <: AbstractPolyhedraOptimizer{T}
     status::MOI.TerminationStatusCode
     solution::Union{AbstractVector{T}, Nothing}
 
+    tol::Union{T,MA.Zero}
+
     function VRepOptimizer{T}(library::Union{Nothing, Library} = nothing) where T
         new(library, LPHRep(_MOIModel{T}()), nothing, nothing,
             MOI.FEASIBILITY_SENSE, nothing, zero(T),
-            MOI.OPTIMIZE_NOT_CALLED, nothing)
+            MOI.OPTIMIZE_NOT_CALLED, nothing,
+            _default_tol(T),
+        )
     end
 end
 
@@ -90,11 +94,11 @@ 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) : _lt(a, b))
+        _better(a, b) = (lpm.objective_sense == MOI.MAX_SENSE ? _gt(a, b; tol) : _lt(a, b; tol))
         bestobjval = zero(T)
         lpm.solution = nothing
         for r in allrays(prob)
-            objval = obj ⋅ r
+            objval = _dot(obj, r)
             if _better(objval, bestobjval)
                 bestobjval = objval
                 lpm.solution = coord(r)
@@ -104,7 +108,7 @@ function MOI.optimize!(lpm::VRepOptimizer{T}) where T
             lpm.status = MOI.DUAL_INFEASIBLE
         else
             for p in points(prob)
-                objval = obj ⋅ p
+                objval = _dot(obj, p)
                 if lpm.solution === nothing || better(objval, bestobjval)
                     bestobjval = objval
                     lpm.solution = p
@@ -137,9 +141,9 @@ function MOI.get(lpm::VRepOptimizer{T}, attr::MOI.ObjectiveValue) where T
         return zero(T)
     end
     if lpm.status == MOI.OPTIMAL
-        return lpm.objective_func ⋅ lpm.solution + lpm.objective_constant
+        return _dot(lpm.objective_func, lpm.solution) + lpm.objective_constant
     elseif lpm.status == MOI.DUAL_INFEASIBLE
-        return lpm.objective_func ⋅ lpm.solution
+        return _dot(lpm.objective_func, lpm.solution)
     else
         error("No objective value available when termination status is $(lpm.status).")
     end
diff --git a/test/comp.jl b/test/comp.jl
index de8b77c4..db65f410 100644
--- a/test/comp.jl
+++ b/test/comp.jl
@@ -1,6 +1,6 @@
 @testset "Comparison" begin
     #@test SymPoint([1, 2]) ≈ SymPoint([-1, -2])
     #@test !(SymPoint([1, 2]) ≈ SymPoint([-2, -1]))
-    @test Polyhedra._lt(1, 1.5)
-    @test !Polyhedra._lt(2, 1.5)
+    @test Polyhedra._lt(1, 1.5; tol = 1e-6)
+    @test !Polyhedra._lt(2, 1.5; tol = 1e-6)
 end
diff --git a/test/issue301.jl b/test/issue301.jl
index c22e7a8f..1eb290fa 100644
--- a/test/issue301.jl
+++ b/test/issue301.jl
@@ -3,8 +3,8 @@ Test for https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/301
 """
 function issue301test(lib::Polyhedra.Library)
     p = polyhedron(vrep([[3, 0], [0, 3], [0, 0], [1, 1]]), lib)
-    removevredundancy!(p; verbose=0, ztol=1e-7)
+    removevredundancy!(p; verbose=0, tol=1e-7)
     @test npoints(p) == 3
-    removehredundancy!(p; verbose=0, ztol=1e-7)
+    removehredundancy!(p; verbose=0, tol=1e-7)
     @test nhalfspaces(p) == 3
 end
diff --git a/test/lp_to_polyhedra.jl b/test/lp_to_polyhedra.jl
index e9e5410e..f0774959 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 mock.objective_func ⋅ mock.solution + mock.objective_constant
+        return _dot(mock.objective_func, mock.solution) + mock.objective_constant
     elseif mock.status == MOI.DUAL_INFEASIBLE
-        return mock.objective_func ⋅ mock.solution
+        return _dot(mock.objective_func, mock.solution)
     else
         error("No objective value available when termination status is $(mock.status).")
     end
diff --git a/test/lphrep.jl b/test/lphrep.jl
index 27572e4d..fbee4dd1 100644
--- a/test/lphrep.jl
+++ b/test/lphrep.jl
@@ -21,7 +21,6 @@ function test_order_lphrep()
     @variable(model, x)
     @variable(model, y <= 1)
     h = hrep(model)
-    @show h
     @test nhalfspaces(h) == 1
     @test first(halfspaces(h)) == HalfSpace([0, 1], 1)
 end
diff --git a/test/redundancy.jl b/test/redundancy.jl
index c4b32dac..4bd4efca 100644
--- a/test/redundancy.jl
+++ b/test/redundancy.jl
@@ -358,13 +358,13 @@ end
         1210.55555992692 -37126.858723322744 -13.95377106842899 -0.0 -0.0 -0.0 -0.0 -0.0
         -1860.676468633397 57155.38024488611 18.21165578673308 -0.0 -0.0 -0.0 -0.0 -0.0
     ]
-    # The last row is almost considered an equality, depends on `ztol`.
+    # The last row is almost considered an equality, depends on `tol`.
     # If it's detected as an equality then the polyhedron is detected as empty.
     # See https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/270
     b = [-20.995, -1.007487377909021, -0.9877340129555137, -0.9705208833525414, -0.9548402631790559, -0.025711479466835108, 1.007487377909021, 0.9877340129555137, 0.9705208833525414, 0.9548402631790559, 0.025711479466835108, 0.8046759919340509, -11751.33121007926, 18148.933864538412]
     lib = DefaultLibrary{Float64}(lp_solver)
     P = polyhedron(hrep(A, b), lib)
-    removevredundancy!(P, ztol=1e-4)
+    removevredundancy!(P, tol=1e-4)
     @test nhyperplanes(P) >= 4
     @test nhalfspaces(P) <= 4
     @test npoints(P) > 0