diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 3f0b2117b..e2b509e45 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -22,6 +22,7 @@ import Base:
     in, zero, one, eps, typemin, typemax, abs, abs2, real, min, max,
     sqrt, exp, log, sin, cos, tan, cot, inv, cbrt, csc, hypot, sec,
     exp2, exp10, log2, log10,
+    mod,
     asin, acos, atan,
     sinh, cosh, tanh, coth, csch, sech, asinh, acosh, atanh, sinpi, cospi,
     union, intersect, isempty,
diff --git a/src/intervals/functions.jl b/src/intervals/functions.jl
index 110bb41ab..aa3ee1c42 100644
--- a/src/intervals/functions.jl
+++ b/src/intervals/functions.jl
@@ -373,3 +373,21 @@ function nthroot(a::Interval{T}, n::Integer) where T
     b = nthroot(bigequiv(a), n)
     return convert(Interval{T}, b)
 end
+
+"""
+Calculate `x::Interval mod y::Real` where y != zero(y) and y is not inteval`.
+"""
+function mod(x::Interval, y::Real)
+    @assert y != zero(y) """mod(x::Interval, y::Real) 
+is currently implemented only for a strictly positive or negative divisor y."""
+    division = x / y
+    fl = floor(division)
+    if !isthin(fl)
+        return y > zero(y) ? Interval(zero(y), y) : Interval(y, zero(y))
+    else
+        return y * (division - fl)
+    end
+end
+
+mod(x::Interval, y::Interval) = throw(ArgumentError("mod not defined for interval as divisor `y`"))
+mod(x::Real, y::Interval) = throw(ArgumentError("mod not defined for interval as divisor `y`"))
diff --git a/test/interval_tests/numeric.jl b/test/interval_tests/numeric.jl
index 3f447a5d6..11f1d0877 100644
--- a/test/interval_tests/numeric.jl
+++ b/test/interval_tests/numeric.jl
@@ -434,3 +434,44 @@ end
     @test nthroot(Interval{BigFloat}(-81, -16), -4) == ∅
     @test nthroot(Interval{BigFloat}(-81, -16), 1) == Interval{BigFloat}(-81, -16)
 end
+
+# approximation used for testing (not to rely on ≈ for intervals)
+# ⪆(x, y) = (x ≈ y) && (x ⊇ y)
+⪆(x::Interval, y::Interval) = x.lo ≈ y.lo && x.hi ≈ y.hi && x ⊇ y
+
+@testset "`mod`" begin
+    r = 0.0625
+    x = r..(1+r)
+    @test mod(x, 1) == mod(x, 1.0) == 0..1
+    @test mod(x, 2) == mod(x, 2.0) ⪆ x
+    @test mod(x, 2.5) ⪆ x
+    @test mod(x, 0.5) == 0..0.5
+    @test mod(x, -1) == mod(x, -1.0) == -1..0
+    @test mod(x, -2) == mod(x, -2.0) ⪆ -2+x
+    @test mod(x, -2.5) ⪆ -2.5+x
+    @test mod(x, -0.5) == -0.5..0
+
+    x = (-1+r) .. -r
+    @test mod(x, 1) == mod(x, 1.0) ⪆ 1+x
+    @test mod(x, 2) == mod(x, 2.0) ⪆ 2+x
+    @test mod(x, 2.5) ⪆ 2.5+x
+    @test mod(x, 0.5) == 0..0.5
+    @test mod(x, -1) == mod(x, -1.0) ⪆ x
+    @test mod(x, -2) == mod(x, -2.0) ⪆ x
+    @test mod(x, -2.5) ⪆ x
+    @test mod(x, -0.5) == -0.5..0
+
+    x = -r .. 1-r
+    @test mod(x, 1) == mod(x, 1.0) == 0..1
+    @test mod(x, 2) == mod(x, 2.0) == 0..2
+    @test mod(x, 2.5) == 0..2.5
+    @test mod(x, 0.5) == 0..0.5
+    @test mod(x, -1) == mod(x, -1.0) == -1..0
+    @test mod(x, -2) == mod(x, -2.0) == -2..0
+    @test mod(x, -2.5) == -2.5..0
+    @test mod(x, -0.5) == -0.5..0
+
+    # TODO - implement mod for two intervals
+    @test_throws ArgumentError mod(1..2, 1.4..1.5)
+    @test_throws ArgumentError mod(1.0, 1.4..1.5)
+end