diff --git a/Project.toml b/Project.toml
index 3dfb05e6..1647abe6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,14 +1,15 @@
 name = "IntervalArithmetic"
 uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-repo = "https://github.com/JuliaIntervals/IntervalArithmetic.jl.git"
 version = "1.0.2"
+repo = "https://github.com/JuliaIntervals/IntervalArithmetic.jl.git"
 
 [deps]
 CRlibm = "96374032-68de-5a5b-8d9e-752f78720389"
+CoreMath = "b7a15901-be09-4a0e-87d2-2e66b0e09b5a"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 OpenBLASConsistentFPCSR_jll = "6cdc7f73-28fd-5e50-80fb-958a8875b1af"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RoundingEmulator = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
 
 [weakdeps]
@@ -32,14 +33,15 @@ IntervalArithmeticSparseArraysExt = "SparseArrays"
 [compat]
 Arblib = "1.3.0"
 CRlibm = "1.0.2"
+CoreMath = "0.1.0"
 DiffRules = "1"
 ForwardDiff = "1"
 IntervalSets = "0.7"
 LinearAlgebra = "1.10"
 MacroTools = "0.5"
 OpenBLASConsistentFPCSR_jll = "0.3.29"
-Random = "1.10"
 Printf = "1.10"
+Random = "1.10"
 RecipesBase = "1"
 RoundingEmulator = "0.2"
 SparseArrays = "1.10"
diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index d12f8141..d837bec2 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -44,7 +44,7 @@ with decorations and up to 6 significant digits.
 """
 module IntervalArithmetic
 
-import RoundingEmulator, CRlibm, Base.MPFR
+import RoundingEmulator, CRlibm, CoreMath, Base.MPFR
 
 using MacroTools: MacroTools, prewalk, postwalk, @capture
 
@@ -102,7 +102,7 @@ function configure_flavor(flavor::Symbol)
 end
 
 function configure_rounding(rounding::Symbol)
-    rounding ∈ (:correct, :none) || return throw(ArgumentError("only the rounding mode `:correct` and `:none` are available"))
+    rounding ∈ (:correct, :ulp, :none) || return throw(ArgumentError("only the rounding mode `:correct`, `:ulp` and `:none` are available"))
     @eval default_rounding() = IntervalRounding{$(QuoteNode(rounding))}()
     return rounding
 end
diff --git a/src/intervals/rounding.jl b/src/intervals/rounding.jl
index 4b8bf295..e84142e2 100644
--- a/src/intervals/rounding.jl
+++ b/src/intervals/rounding.jl
@@ -4,8 +4,8 @@
 Interval rounding type.
 
 Available rounding types:
-- `:correct`: rounding via [RoundingEmulator.jl](https://github.com/JuliaIntervals/RoundingEmulator.jl)
-and [CRlibm.jl](https://github.com/JuliaInterval/IntervalArithmetic.jl).
+- `:correct`: rounding via [RoundingEmulator.jl](https://github.com/JuliaIntervals/RoundingEmulator.jl) and [CRlibm.jl](https://github.com/JuliaInterval/IntervalArithmetic.jl); fallback to MPFR.
+- `:ulp`: rounding via [CoreMath.jl](https://github.com/JuliaMath/CoreMath.jl) (default rounding to nearest) with `prevfloat` and `nextfloat`; fallback to MPFR.
 - `:none`: no rounding (non-rigorous numerics).
 """
 struct IntervalRounding{T} end
@@ -15,22 +15,21 @@ struct IntervalRounding{T} end
 _fround(f::Function, x, y, r) = _fround(f, default_rounding(), x, y, r)
 _fround(f::Function, x, r)    = _fround(f, default_rounding(), x, r)
 
-#
+# +, -, *, /, inv, sqrt
 
 for (f, fname) ∈ ((:+, :add), (:-, :sub), (:*, :mul), (:/, :div))
+
+    mpfr_f = Symbol(:mpfr_, fname)
+
     @eval begin
         _fround(::typeof($f), ::IntervalRounding, x::T, y::T, ::RoundingMode) where {T<:Rational} = $f(x, y) # exact operation
 
-        _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Down}) where {T<:Union{Float32,Float64}} =
-            RoundingEmulator.$(Symbol(fname, :_down))(x, y)
-        _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} =
-            RoundingEmulator.$(Symbol(fname, :_up))(x, y)
-        function _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
+        function _fround(::typeof($f), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
             prec = max(precision(x), precision(y))
             bigx = BigFloat(x; precision = prec)
             bigy = BigFloat(y; precision = prec)
             bigz = BigFloat(; precision = prec)
-            @ccall Base.MPFR.libmpfr.$(Symbol(:mpfr_, fname))(
+            @ccall Base.MPFR.libmpfr.$mpfr_f(
                 bigz::Ref{BigFloat},
                 bigx::Ref{BigFloat},
                 bigy::Ref{BigFloat},
@@ -39,64 +38,29 @@ for (f, fname) ∈ ((:+, :add), (:-, :sub), (:*, :mul), (:/, :div))
             return bigz
         end
 
-        _fround(::typeof($f), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = $f(x, y)
-    end
-end
-
-#
+        _fround(::typeof($f), ir::IntervalRounding{:correct}, x::Float16, y::Float16, r::RoundingMode) =
+            Float16(_fround($f, ir, Float64(x), Float64(y), r), r)
+        _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Down}) where {T<:Union{Float32,Float64}} =
+            RoundingEmulator.$(Symbol(fname, :_down))(x, y)
+        _fround(::typeof($f), ::IntervalRounding{:correct}, x::T, y::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} =
+            RoundingEmulator.$(Symbol(fname, :_up))(x, y)
 
-_fround(::typeof(^), ::IntervalRounding, x::T, y::T, ::RoundingMode) where {T<:Rational} = ^(x, y) # exact operation
-_fround(::typeof(^), ::IntervalRounding, x::Rational, n::Integer, ::RoundingMode) = ^(x, n) # exact operation
+        _fround(::typeof($f), ir::IntervalRounding{:ulp}, x::T, y::T, r::RoundingMode) where {T<:Union{Float16,Float32}} =
+            T(_fround($f, ir, Float64(x), Float64(y), r), r)
+        _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Float64, y::Float64, ::RoundingMode{:Down}) = prevfloat($f(x, y))
+        _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Float64, y::Float64, ::RoundingMode{:Up}) = nextfloat($f(x, y))
 
-_fround(::typeof(^), ir::IntervalRounding, x::AbstractFloat, n::Integer, r::RoundingMode) =
-    _fround(^, ir, promote(x, n)..., r)
-
-function _fround(::typeof(^), ::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
-    prec = max(precision(x), precision(y))
-    bigx = BigFloat(x; precision = prec)
-    bigy = BigFloat(y; precision = prec)
-    bigz = BigFloat(; precision = prec)
-    @ccall Base.MPFR.libmpfr.mpfr_pow(
-        bigz::Ref{BigFloat},
-        bigx::Ref{BigFloat},
-        bigy::Ref{BigFloat},
-        r::Base.MPFR.MPFRRoundingMode
-    )::Int32
-    return bigz
+        _fround(::typeof($f), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = $f(x, y)
+    end
 end
 
-_fround(::typeof(^), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = ^(x, y)
-
-#
-
 _fround(::typeof(inv), ::IntervalRounding, x::Rational, ::RoundingMode) = inv(x) # exact operation
 
-_fround(::typeof(inv), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) =
-    RoundingEmulator.div_down(one(x), x)
-_fround(::typeof(inv), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) =
-    RoundingEmulator.div_up(one(x), x)
-function _fround(::typeof(inv), ::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode)
-    prec = precision(x)
-    bigx = BigFloat(x; precision = prec)
-    bigz = BigFloat(; precision = prec)
-    @ccall Base.MPFR.libmpfr.mpfr_div(
-        bigz::Ref{BigFloat},
-        one(bigx)::Ref{BigFloat},
-        bigx::Ref{BigFloat},
-        r::Base.MPFR.MPFRRoundingMode
-    )::Int32
-    return bigz
-end
-
-_fround(::typeof(inv), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = inv(x)
+_fround(::typeof(inv), ::IntervalRounding, x::AbstractFloat, r::RoundingMode) = _fround(/, one(x), x, r)
 
-#
+_fround(::typeof(inv), ::IntervalRounding{:none}, x::T, ::RoundingMode) where {T<:AbstractFloat} = inv(x)
 
-_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Down}) =
-    RoundingEmulator.sqrt_down(x)
-_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, ::RoundingMode{:Up}) =
-    RoundingEmulator.sqrt_up(x)
-function _fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode)
+function _fround(::typeof(sqrt), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::T, r::RoundingMode) where {T<:AbstractFloat}
     prec = precision(x)
     bigx = BigFloat(x; precision = prec)
     bigz = BigFloat(; precision = prec)
@@ -108,68 +72,33 @@ function _fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::AbstractFloat,
     return bigz
 end
 
-_fround(::typeof(sqrt), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = sqrt(x)
-
-#
-
-function rootn end # defined in arithmetic/power.jl
-
-function _fround(::typeof(rootn), ::IntervalRounding{:correct}, x::AbstractFloat, n::Integer, r::RoundingMode)
-    prec = precision(x)
-    bigx = BigFloat(x; precision = prec)
-    bigz = BigFloat(; precision = prec)
-    @ccall Base.MPFR.libmpfr.mpfr_rootn_ui(
-        bigz::Ref{BigFloat},
-        bigx::Ref{BigFloat},
-        n::Culong,
-        r::Base.MPFR.MPFRRoundingMode
-    )::Int32
-    return bigz
-end
-
-_fround(::typeof(rootn), ::IntervalRounding{:none}, x::AbstractFloat, n::Integer, ::RoundingMode) = x^(1//n)
-
-#
+_fround(::typeof(sqrt), ir::IntervalRounding{:correct}, x::Float16, r::RoundingMode) =
+    Float16(_fround(sqrt, ir, Float64(x), r), r)
+_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::T, ::RoundingMode{:Down}) where {T<:Union{Float32,Float64}} =
+    RoundingEmulator.sqrt_down(x)
+_fround(::typeof(sqrt), ::IntervalRounding{:correct}, x::T, ::RoundingMode{:Up}) where {T<:Union{Float32,Float64}} =
+    RoundingEmulator.sqrt_up(x)
 
-function _fround(::typeof(atan), ::IntervalRounding{:correct}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
-    prec = max(precision(x), precision(y))
-    bigx = BigFloat(x; precision = prec)
-    bigy = BigFloat(y; precision = prec)
-    bigz = BigFloat(; precision = prec)
-    @ccall Base.MPFR.libmpfr.mpfr_atan2(
-        bigz::Ref{BigFloat},
-        bigx::Ref{BigFloat},
-        bigy::Ref{BigFloat},
-        r::Base.MPFR.MPFRRoundingMode
-    )::Int32
-    return bigz
-end
+_fround(::typeof(sqrt), ir::IntervalRounding{:ulp}, x::T, r::RoundingMode) where {T<:Union{Float16,Float32}} =
+    T(_fround(sqrt, ir, Float64(x), r), r)
+_fround(::typeof(sqrt), ::IntervalRounding{:ulp}, x::Float64, ::RoundingMode{:Down}) = prevfloat(sqrt(x))
+_fround(::typeof(sqrt), ::IntervalRounding{:ulp}, x::Float64, ::RoundingMode{:Up}) = nextfloat(sqrt(x))
 
-_fround(::typeof(atan), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = atan(x, y)
+_fround(::typeof(sqrt), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = sqrt(x)
 
-#
+# 1-argument functions
 
-for f ∈ (:cbrt, :exp2, :exp10, :cot, :sec, :csc, :tanh, :coth, :sech, :csch, :asinh, :acosh, :atanh)
-    @eval begin
-        function _fround(::typeof($f), ::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode)
-            prec = precision(x)
-            bigx = BigFloat(x; precision = prec)
-            bigz = BigFloat(; precision = prec)
-            @ccall Base.MPFR.libmpfr.$(Symbol(:mpfr_, f))(
-                bigz::Ref{BigFloat},
-                bigx::Ref{BigFloat},
-                r::Base.MPFR.MPFRRoundingMode
-            )::Int32
-            return bigz
-        end
+for f ∈ (:cbrt, :exp, :exp2, :exp10, :expm1, # exponential
+         :log, :log2, :log10, :log1p, # logarithmic
+         :sin, :sinpi, :cos, :cospi, :tan, :cot, :sec, :csc, :asin, :acos, :atan, :acot, # trigonometric
+         :sinh, :tanh, :asinh, :cosh, :coth, :sech, :csch, :acosh, :atanh, :acoth) # hyperbolic
 
-        _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x)
-    end
-end
+    mpfr_f = Symbol(:mpfr_, f)
+    coremath_f = Symbol(:cr_, f)
 
-for (f, g) ∈ ((:acot, :atan), (:acoth, :atanh))
-    @eval begin
-        function _fround(::typeof($f), ir::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode{:Down})
+    if f ∈ (:acot, :acoth)
+        g = f === :acot ? :atan : :atanh
+        @eval function _fround(::typeof($f), ir::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, r::RoundingMode{:Down})
             prec = precision(x)
             bigx = BigFloat(x; precision = prec + 10)
             bigz = BigFloat(; precision = prec + 10)
@@ -182,7 +111,7 @@ for (f, g) ∈ ((:acot, :atan), (:acoth, :atanh))
             bigw = _fround($g, ir, bigz, r)
             return BigFloat(bigw, r; precision = prec)
         end
-        function _fround(::typeof($f), ir::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode{:Up})
+        @eval function _fround(::typeof($f), ir::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, r::RoundingMode{:Up})
             prec = precision(x)
             bigx = BigFloat(x; precision = prec + 32)
             bigz = BigFloat(; precision = prec + 32)
@@ -195,23 +124,82 @@ for (f, g) ∈ ((:acot, :atan), (:acoth, :atanh))
             bigw = _fround($g, ir, bigz, r)
             return BigFloat(bigw, r; precision = prec)
         end
+    else
+        @eval function _fround(::typeof($f), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, r::RoundingMode)
+            prec = precision(x)
+            bigx = BigFloat(x; precision = prec)
+            bigz = BigFloat(; precision = prec)
+            @ccall Base.MPFR.libmpfr.$mpfr_f(
+                bigz::Ref{BigFloat},
+                bigx::Ref{BigFloat},
+                r::Base.MPFR.MPFRRoundingMode
+            )::Int32
+            return bigz
+        end
+    end
 
-        _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, ::RoundingMode) = $f(x)
+    if f ∈ CRlibm.functions
+        @eval _fround(::typeof($f), ::IntervalRounding{:correct}, x::Float16, r::RoundingMode) = Float16(_fround($f, Float64(x), r), r)
+        @eval _fround(::typeof($f), ::IntervalRounding{:correct}, x::Union{Float32,Float64}, r::RoundingMode) = CRlibm.$f(x, r)
     end
+
+    @eval _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Union{Float16,Float32,Float64}, ::RoundingMode{:Down}) = prevfloat(CoreMath.$coremath_f(x))
+    @eval _fround(::typeof($f), ::IntervalRounding{:ulp}, x::Union{Float16,Float32,Float64}, ::RoundingMode{:Up}) = nextfloat(CoreMath.$coremath_f(x))
+
+    @eval _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, r::RoundingMode) = $f(x)
 end
 
-# CRlibm functions
+# 2-argument functions
+
+_fround(::typeof(^), ::IntervalRounding, x::Rational, n::Integer, ::RoundingMode) = ^(x, n) # exact operation
+
+_fround(::typeof(^), ir::IntervalRounding, x::AbstractFloat, n::Integer, r::RoundingMode) =
+    _fround(^, ir, promote(x, n)..., r)
+
+for (f, fname) ∈ ((:^, :pow), (:atan, :atan2))
 
-for f ∈ CRlibm.functions
-    if f !== :atanpi # currently not defined in IntervalArithmetic
-        @eval begin
-            _fround(::typeof($f), ::IntervalRounding{:correct}, x::AbstractFloat, r::RoundingMode) = CRlibm.$f(x, r)
+    mpfr_f = Symbol(:mpfr_, fname)
+    coremath_f = Symbol(:cr_, fname)
 
-            _fround(::typeof($f), ::IntervalRounding{:none}, x::AbstractFloat, r::RoundingMode) = $f(x)
+    @eval begin
+        function _fround(::typeof($f), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::T, y::T, r::RoundingMode) where {T<:AbstractFloat}
+            prec = max(precision(x), precision(y))
+            bigx = BigFloat(x; precision = prec)
+            bigy = BigFloat(y; precision = prec)
+            bigz = BigFloat(; precision = prec)
+            @ccall Base.MPFR.libmpfr.$mpfr_f(
+                bigz::Ref{BigFloat},
+                bigx::Ref{BigFloat},
+                bigy::Ref{BigFloat},
+                r::Base.MPFR.MPFRRoundingMode
+            )::Int32
+            return bigz
         end
+
+        _fround(::typeof($f), ::IntervalRounding{:ulp}, x::T, y::T, ::RoundingMode{:Down}) where {T<:Union{Float16,Float32,Float64}} = prevfloat(CoreMath.$coremath_f(x, y))
+        _fround(::typeof($f), ::IntervalRounding{:ulp}, x::T, y::T, ::RoundingMode{:Up}) where {T<:Union{Float16,Float32,Float64}} = nextfloat(CoreMath.$coremath_f(x, y))
+
+        _fround(::typeof($f), ::IntervalRounding{:none}, x::T, y::T, ::RoundingMode) where {T<:AbstractFloat} = $f(x, y)
     end
 end
 
+function rootn end # defined in arithmetic/power.jl
+
+function _fround(::typeof(rootn), ::Union{IntervalRounding{:correct},IntervalRounding{:ulp}}, x::AbstractFloat, n::Integer, r::RoundingMode)
+    prec = precision(x)
+    bigx = BigFloat(x; precision = prec)
+    bigz = BigFloat(; precision = prec)
+    @ccall Base.MPFR.libmpfr.mpfr_rootn_ui(
+        bigz::Ref{BigFloat},
+        bigx::Ref{BigFloat},
+        n::Culong,
+        r::Base.MPFR.MPFRRoundingMode
+    )::Int32
+    return bigz
+end
+
+_fround(::typeof(rootn), ::IntervalRounding{:none}, x::AbstractFloat, n::Integer, ::RoundingMode) = x^(1//n)
+