diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index 910ca1b48..3639e8924 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -27,6 +27,7 @@ import Base:
     in, zero, one, eps, typemin, typemax, abs, abs2, real, min, max,
     sqrt, exp, log, sin, cos, tan, inv,
     exp2, exp10, log2, log10,
+    mod,
     asin, acos, atan, atan2,
     sinh, cosh, tanh, asinh, acosh, atanh,
     union, intersect, isempty,
@@ -61,7 +62,7 @@ export
     RoundTiesToEven, RoundTiesToAway,
     cancelminus, cancelplus, isunbounded,
     .., @I_str, ±,
-    pow, extended_div,
+    pow, extended_div, extended_mod,
     setformat, @format
 
 export
diff --git a/src/intervals/functions.jl b/src/intervals/functions.jl
index e810afc0a..3a40fa6fb 100644
--- a/src/intervals/functions.jl
+++ b/src/intervals/functions.jl
@@ -266,3 +266,70 @@ for f in (:log, :log2, :log10, :log1p)
 
         end
 end
+
+# mod
+
+"""
+Calculate `X mod a` where `X` is an interval and `a` is a positive, atomic interval.
+"""
+function mod(X::Interval, a::Interval)
+    division = X / a
+    fl = floor(division)
+
+    if fl.lo < fl.hi
+        return 0..(a.hi)
+    end
+
+    return (division - fl) * a
+end
+
+mod(X::Interval, a::Real) = mod(X, interval(a))
+#
+# function mod(a::Interval, y::T) where {T<:Real}
+#     yy = abs(y)
+#     fld_lo = floor(a.lo/yy)
+#     fld_hi = floor(a.hi/yy)
+#     z = zero(fld_lo)
+#
+#     if fld_lo != fld_hi
+#         # `a` includes a discontinuity of `mod`
+#         if y > 0
+#             return interval(z, y)
+#         else
+#             return interval(y, z)
+#         end
+#     else
+#         # no discontinuity crossed within `a`
+#         return interval(mod(a.lo, y), mod(a.hi, y))
+#     end
+# end
+
+
+function extended_mod(a::Interval, y::T) where {T<:Real}
+    yy = abs(y)
+    fld_lo = floor(a.lo/yy)
+    fld_hi = floor(a.hi/yy)
+    z = zero(fld_lo)
+    S = typeof( z )
+    ee = emptyinterval(S)
+
+    if fld_lo + 1 == fld_hi
+        # `a` includes one discontinuity
+        if y > 0
+            return interval(mod(a.lo, y), y), interval(z, mod(a.hi, y)), ee
+        else
+            return interval(mod(a.lo, y), z), interval(y, mod(a.hi, y)), ee
+        end
+    elseif fld_lo +1 < fld_hi
+        # `a` includes more discontinuities
+        if y > 0
+            return interval(mod(a.lo, y), y), interval(z, y), interval(z, mod(a.hi, y))
+        else
+            return interval(mod(a.lo, y), z), interval(y, z), interval(y, mod(a.hi, y))
+        end
+    else
+        # no discontinuity crossed within `a`
+        return interval(mod(a.lo, y), mod(a.hi, y)), ee, ee
+    end
+
+end
diff --git a/test/interval_tests/numeric.jl b/test/interval_tests/numeric.jl
index 82c321e88..18c39ff12 100644
--- a/test/interval_tests/numeric.jl
+++ b/test/interval_tests/numeric.jl
@@ -284,3 +284,85 @@ end
     end
     =#
 end
+
+@testset "`mod`" begin
+    r = 0.0625
+    a = r..(1+r)
+    @test mod(a, 1) == mod(a, 1.0) == 0..1
+    @test mod(a, 2) == mod(a, 2.0) == a
+    @test mod(a, 2.5) == a
+    @test mod(a, 0.5) == 0..0.5
+    @test mod(a, -1) == mod(a, -1.0) == -1..0
+    @test mod(a, -2) == mod(a, -2.0) == -2+a
+    @test mod(a, -2.5) == -2.5+a
+    @test mod(a, -0.5) == -0.5..0
+
+    a = (-1+r) .. -r
+    @test mod(a, 1) == mod(a, 1.0) == 1+a
+    @test mod(a, 2) == mod(a, 2.0) == 2+a
+    @test mod(a, 2.5) == 2.5+a
+    @test mod(a, 0.5) == 0..0.5
+    @test mod(a, -1) == mod(a, -1.0) == a
+    @test mod(a, -2) == mod(a, -2.0) == a
+    @test mod(a, -2.5) == a
+    @test mod(a, -0.5) == -0.5..0
+
+    a = -r .. 1-r
+    @test mod(a, 1) == mod(a, 1.0) == 0..1
+    @test mod(a, 2) == mod(a, 2.0) == 0..2
+    @test mod(a, 2.5) == 0..2.5
+    @test mod(a, 0.5) == 0..0.5
+    @test mod(a, -1) == mod(a, -1.0) == -1..0
+    @test mod(a, -2) == mod(a, -2.0) == -2..0
+    @test mod(a, -2.5) == -2.5..0
+    @test mod(a, -0.5) == -0.5..0
+
+    a = pi_interval(Float64)
+    @test mod(a, pi) == interval(0.0, a.hi-pi)
+    @test mod(2a, pi) == interval(0.0, 2*(a.hi-pi))
+    @test mod(1.5a, pi) == interval(0.5a.lo, 1.5a.hi-pi)
+    @test mod(interval(0.25,3pi), pi) == interval(0.0, a.lo)
+end
+
+@testset "`extended_mod`" begin
+    r = 0.0625
+    a = r..(1+r)
+    ee = emptyinterval(Float64)
+    @test extended_mod(a, 1) == (r..1, 0..r, ee)
+    @test extended_mod(a, 2) == (a, ee, ee)
+    @test extended_mod(a, 2.5) == (a, ee, ee)
+    @test extended_mod(a, 0.5) == (r..0.5, 0..0.5, 0..r)
+    @test extended_mod(a, -1) == ((-1+r)..0, -1..(-1+r), ee)
+    @test extended_mod(a, -2) == ((-2+r)..(-1+r), ee, ee)
+    @test extended_mod(a, -2.5) == ((-2.5+r)..(-1.5+r), ee, ee)
+    @test extended_mod(a, -0.5) == ((-0.5+r)..0, -0.5..0, -0.5..(-0.5+r))
+
+    a = (-1+r) .. -r
+    @test extended_mod(a, 1) == (1+a, ee, ee)
+    @test extended_mod(a, 2) == (2+a, ee, ee)
+    @test extended_mod(a, 2.5) == (2.5+a, ee, ee)
+    @test extended_mod(a, 0.5) == (r..0.5, 0..(0.5-r), ee)
+    @test extended_mod(a, -1) == ((-1+r) .. -r, ee, ee)
+    @test extended_mod(a, -2) == ((-1+r) .. -r, ee, ee)
+    @test extended_mod(a, -2.5) == ((-1+r) .. -r, ee, ee)
+    @test extended_mod(a, -0.5) == ((-0.5+r)..0, -0.5 .. -r, ee)
+
+    a = -r .. 1-r
+    @test extended_mod(a, 1) == ((1-r)..1, 0..(1-r), ee)
+    @test extended_mod(a, 2) == ((2-r)..2, 0..(1-r), ee)
+    @test extended_mod(a, 2.5) == ((2.5-r)..2.5, 0..(1-r), ee)
+    @test extended_mod(a, 0.5) == ((0.5-r)..0.5, 0..0.5, 0..(0.5-r))
+    @test extended_mod(a, -1) == (-r..0, -1..(-r), ee)
+    @test extended_mod(a, -2) == (-r..0, -2..(-1-r), ee)
+    @test extended_mod(a, -2.5) == (-r..0, (-2.5)..(-1.5-r), ee)
+    @test extended_mod(a, -0.5) == (-r..0, -0.5..0, -0.5..(-r))
+
+    a = pi_interval(Float64)
+    @test extended_mod(a, pi) == (mod(a, pi), ee, ee)
+    @test extended_mod(2a, pi) == (interval(0.0, 2*(a.hi-pi)), ee, ee)
+    @test extended_mod(1.5a, pi) == (interval(0.5a.lo, 1.5a.hi-pi), ee, ee)
+    @test extended_mod(interval(0.25,3pi), pi) == (interval(0.25, a.lo),
+        interval(0.0, a.lo), interval(0.0))
+    @test extended_mod(@interval(0.25,3pi), pi) == (interval(0.25, a.lo),
+        interval(0.0, a.lo), mod(@interval(3pi), pi))
+end