From f886a50ffc73dd89c6b702f2e00256fdfd9e6be4 Mon Sep 17 00:00:00 2001
From: Luis Benet <benet@fis.unam.mx>
Date: Fri, 4 May 2018 10:43:32 -0500
Subject: [PATCH] Add extended_mod(a, y), with tests

---
 src/IntervalArithmetic.jl      |  2 +-
 src/intervals/functions.jl     | 30 ++++++++++++++++++
 test/interval_tests/numeric.jl | 56 +++++++++++++++++++++++++++++-----
 3 files changed, 79 insertions(+), 9 deletions(-)

diff --git a/src/IntervalArithmetic.jl b/src/IntervalArithmetic.jl
index c2df13d15..bee24876d 100644
--- a/src/IntervalArithmetic.jl
+++ b/src/IntervalArithmetic.jl
@@ -45,7 +45,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 08eaf6ee0..534c68746 100644
--- a/src/intervals/functions.jl
+++ b/src/intervals/functions.jl
@@ -286,3 +286,33 @@ function mod(a::Interval, y::T) where {T<:Real}
         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 = fld(a.lo, yy)
+    fld_hi = fld(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 ae5806089..eb652d891 100644
--- a/test/interval_tests/numeric.jl
+++ b/test/interval_tests/numeric.jl
@@ -282,26 +282,66 @@ end
 @testset "`mod`" begin
     r = 0.0625
     a = r..(1+r)
-    @test mod(a, 1) == mod(a, 1.0) == interval(0,1)
+    @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, -1) == mod(a, -1.0) == interval(-1,0)
+    @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) == interval(0,1)
-    @test mod(a, 2) == mod(a, 2.0) == interval(0,2)
-    @test mod(a, 2.5) == interval(0,2.5)
-    @test mod(a, -1) == mod(a, -1.0) == interval(-1,0)
-    @test mod(a, -2) == mod(a, -2.0) == interval(-2,0)
-    @test mod(a, -2.5) == interval(-2.5,0)
+    @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
+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))
 end