diff --git a/Project.toml b/Project.toml
index d94a5953..9693f0af 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Tulip"
 uuid = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"
-authors = ["Mathieu Tanneau <mathieu.tanneau@gmail.com>"]
 version = "0.9.8"
+authors = ["Mathieu Tanneau <mathieu.tanneau@gmail.com>"]
 
 [deps]
 CodecBzip2 = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
diff --git a/src/IPM/ipmdata.jl b/src/IPM/ipmdata.jl
index 7222ada0..04031e27 100644
--- a/src/IPM/ipmdata.jl
+++ b/src/IPM/ipmdata.jl
@@ -95,7 +95,7 @@ function IPMData(pb::ProblemData{T}, mfact::Factory) where{T}
             push!(uslack, T(Inf))
             b[i] = ub
 
-        elseif isfinite(lb) && ub == Inf
+        elseif isfinite(lb) && T(Inf) == ub
             # a'x >= b --> a'x - s = b
             push!(sind, i)
             push!(sval, -one(T))
diff --git a/src/Presolve/Presolve.jl b/src/Presolve/Presolve.jl
index 3778a932..1fa7442c 100644
--- a/src/Presolve/Presolve.jl
+++ b/src/Presolve/Presolve.jl
@@ -173,6 +173,31 @@ mutable struct PresolveData{T}
     end
 end
 
+const IEEEFloat = Union{Float16, Float32, Float64, BigFloat}
+
+infminus(x::T, y::T) where {T<:IEEEFloat} = x - y
+function infminus(x::T, y::T) where {T<:AbstractFloat}
+    isfinite(x) && isfinite(y) && return x - y
+    isfinite(x) && !isfinite(y) && return -y
+    !isfinite(x) && isfinite(y) && return x
+    !isfinite(x) && !isfinite(y) && return signbit(x) != signbit(y) ? x : T(NaN)
+end
+
+inftimes(x::T, y::T) where {T<:IEEEFloat} = x * y
+function inftimes(x::T, y::T) where {T<:AbstractFloat}
+    isfinite(x) && isfinite(y) && return x * y
+    isfinite(x) && !isfinite(y) && return signbit(x) ? -y : y
+    !isfinite(x) && isfinite(y) && return signbit(y) ? -x : x
+    !isfinite(x) && !isfinite(y) && return signbit(x) == signbit(y) ? abs(x) : -abs(x)
+end
+
+infdiv(x::T, y::T) where {T<:IEEEFloat} = x / y
+function infdiv(x::T, y::T) where {T<:AbstractFloat}
+    isfinite(x) && !isfinite(y) && return zero(T)
+    !isfinite(x) && isfinite(y) && return signbit(x) == signbit(y) ? abs(x) : -abs(x)
+    return x / y
+end
+
 # Extract pre-solved problem data, to be passed to the IPM solver
 function extract_reduced_problem!(ps::PresolveData{T}) where{T}
 
@@ -279,8 +304,8 @@ function extract_reduced_problem!(ps::PresolveData{T}) where{T}
             row.nzval[k] /= (rscale[i] * cscale[j])
         end
         # Scale row bounds
-        pb.lcon[i] /= rscale[i]
-        pb.ucon[i] /= rscale[i]
+        pb.lcon[i] = infdiv(pb.lcon[i], rscale[i])
+        pb.ucon[i] = infdiv(pb.ucon[i], rscale[i])
     end
     # Columns
     for (j, col) in enumerate(pb.acols)
diff --git a/src/Presolve/free_column_singleton.jl b/src/Presolve/free_column_singleton.jl
index 21602b8b..180a9eca 100644
--- a/src/Presolve/free_column_singleton.jl
+++ b/src/Presolve/free_column_singleton.jl
@@ -47,30 +47,30 @@ function remove_free_column_singleton!(ps::PresolveData{T}, j::Int) where{T}
                 (ps.colflag[k] && k != j) || continue
                 # Update bounds
                 if aik > 0
-                    l_ -= aik * ps.ucol[k]
-                    u_ -= aik * ps.lcol[k]
+                    l_ = infminus(l_, inftimes(aik, ps.ucol[k]))
+                    u_ = infminus(u_, inftimes(aik, ps.lcol[k]))
                 else
-                    l_ -= aik * ps.lcol[k]
-                    u_ -= aik * ps.ucol[k]
+                    l_ = infminus(l_, inftimes(aik, ps.lcol[k]))
+                    u_ = infminus(u_, inftimes(aik, ps.ucol[k]))
                 end
             end
-            l_ /= aij
-            u_ /= aij
+            l_ = infdiv(l_, aij)
+            u_ = infdiv(u_, aij)
         else
             l_, u_ = ur, lr
             for (k, aik) in zip(row.nzind, row.nzval)
                 (ps.colflag[k] && k != j) || continue
                 # Update bounds
                 if aik > 0
-                    l_ -= aik * ps.lcol[k]
-                    u_ -= aik * ps.ucol[k]
+                    l_ = infminus(l_, inftimes(aik, ps.lcol[k]))
+                    u_ = infminus(u_, inftimes(aik, ps.ucol[k]))
                 else
-                    l_ -= aik * ps.ucol[k]
-                    u_ -= aik * ps.lcol[k]
+                    l_ = infminus(l_, inftimes(aik, ps.ucol[k]))
+                    u_ = infminus(u_, inftimes(aik, ps.lcol[k]))
                 end
             end
-            l_ /= aij
-            u_ /= aij
+            l_ = infdiv(l_, aij)
+            u_ = infdiv(u_, aij)
         end
         @debug """Column singleton $j
             Original bounds: [$l, $u]
diff --git a/test/Project.toml b/test/Project.toml
index b9171ab3..71483068 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,6 +2,7 @@
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -9,3 +10,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 Krylov = "0.10"
 MathOptInterface= "1"
+MultiFloats = "3"
diff --git a/test/runtests.jl b/test/runtests.jl
index ff242181..bdf147f7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,11 +2,12 @@ using LinearAlgebra
 using SparseArrays
 using Test
 using TOML
+using MultiFloats
 
 using Tulip
 TLP = Tulip
 
-const TvTYPES = [Float32, Float64, BigFloat]
+const TvTYPES = (Float32, Float64, Float64x2, BigFloat)
 
 # Check That Tulip.version() matches what's in the Project.toml
 tlp_ver = Tulip.version()