diff --git a/REQUIRE b/REQUIRE
index efc5e6d..f81d81e 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -1,5 +1,6 @@
 julia 0.5
 IntervalArithmetic 0.9
 IntervalRootFinding 0.1
+IntervalContractors 0.1
 MacroTools 0.3
 
diff --git a/src/IntervalConstraintProgramming.jl b/src/IntervalConstraintProgramming.jl
index 936e75a..1c45d5b 100644
--- a/src/IntervalConstraintProgramming.jl
+++ b/src/IntervalConstraintProgramming.jl
@@ -2,7 +2,9 @@ __precompile__()
 
 module IntervalConstraintProgramming
 
-using IntervalArithmetic, IntervalRootFinding
+using   IntervalArithmetic,
+        IntervalRootFinding,
+        IntervalContractors
 
 using MacroTools
 
@@ -20,8 +22,8 @@ export
     Vol,
     show_code
 
+const reverse_operations = IntervalContractors.reverse_operations
 
-include("reverse_mode.jl")
 include("ast.jl")
 include("code_generation.jl")
 include("contractor.jl")
diff --git a/src/ast.jl b/src/ast.jl
index 13bb2b8..1c842a3 100644
--- a/src/ast.jl
+++ b/src/ast.jl
@@ -309,7 +309,7 @@ function process_call!(flatAST::FlatAST, ex, new_var=nothing)
 
     #@show op
 
-    if op ∈ keys(rev_ops)  # standard operator
+    if op ∈ keys(reverse_operations)  # standard operator
         if new_var == nothing
             new_var = make_symbol()
         end
diff --git a/src/code_generation.jl b/src/code_generation.jl
index c08c86e..e4a9474 100644
--- a/src/code_generation.jl
+++ b/src/code_generation.jl
@@ -53,7 +53,7 @@ function emit_backward_code(a::Assignment)
     args = isa(a.args, Vector) ? a.args : [a.args]
 
     return_args = [a.lhs, args...]
-    rev_op = rev_ops[a.op]  # find reverse operation
+    rev_op = reverse_operations[a.op]  # find reverse operation
 
     if rev_op == :()   # empty
         args = make_tuple(args)
diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
deleted file mode 100644
index 185b109..0000000
--- a/src/reverse_mode.jl
+++ /dev/null
@@ -1,130 +0,0 @@
-export plus_rev, minus_rev, mul_rev,
-    power_rev, sqrt_rev, sqr_rev # export due to quoting issue
-
-# export sqr
-# sqr(x) = x^2
-
-const rev_ops = Dict(
-                    :+     => :plus_rev,
-                    :-     => :minus_rev,
-                    :*     => :mul_rev,
-                    :^     => :power_rev,
-                    :sqrt  => :sqrt_rev,
-                    :sqr   => :sqr_rev,
-                    :()    => :()
-                    )
-
-
-function plus_rev(a::Interval, b::Interval, c::Interval)  # a = b + c
-    # a = a ∩ (b + c)  # add this line for plus contractor (as opposed to reverse function)
-    b = b ∩ (a - c)
-    c = c ∩ (a - b)
-
-    return a, b, c
-end
-
-plus_rev(a,b,c) = plus_rev(promote(a,b,c)...)
-
-function minus_rev(a::Interval, b::Interval, c::Interval)  # a = b - c
-    # a = a ∩ (b - c)
-    b = b ∩ (a + c)
-    c = c ∩ (b - a)
-
-    return a, b, c
-end
-
-minus_rev(a,b,c) = minus_rev(promote(a,b,c)...)
-
-minus_rev(a::Interval, b::Interval) = (b = -a; return (a, b))     # a = -b
-
-function mul_rev(a::Interval, b::Interval, c::Interval)  # a = b * c
-    # a = a ∩ (b * c)
-    b = b ∩ (a / c)
-    c = c ∩ (a / b)
-
-    return a, b, c
-end
-
-mul_rev(a,b,c) = mul_rev(promote(a,b,c)...)
-
-
-function power_rev(a::Interval, b::Interval, c::Integer)  # a = b^c,  log(a) = c.log(b),  b = a^(1/c)
-
-    if c == 2  # a = b^2
-        b1 = b ∩ √a
-        b2 = b ∩ (-√a)
-
-        b = hull(b1, b2)
-
-    elseif iseven(c)
-        b1 = b ∩ ( a^(inv(c) ))
-        b2 = b ∩ ( -( a^(inv(c)) ) )
-
-        b = hull(b1, b2)
-
-    elseif isodd(c)
-        b1 = b ∩ ( (a ∩ (0..∞)) ^(inv(c) ))   # positive part
-        b2 = b ∩ (- ( (-(a ∩ (-∞..0)))^(inv(c)) ) )  # negative part
-
-        b = hull(b1, b2)
-    end
-
-    return (a, b, c)
-end
-
-
-function power_rev(a::Interval, b::Interval, c::Interval)  # a = b^c
-
-    # log(a) = c.log(b),  b = a^(1/c)
-
-    b = b ∩ ( a^(inv(c) ))
-    c = c ∩ (log(a) / log(b))
-
-    return a, b, c
-end
-
-power_rev(a, b, c) = power_rev(promote(a, b, c)...)
-
-
-function sqrt_rev(a::Interval, b::Interval)  # a = sqrt(b)
-    # a1 = a ∩ √b
-    # a2 = a ∩ (-(√b))
-    # a = hull(a1, a2)
-
-    b = b ∩ (a^2)
-
-    return a, b
-end
-
-sqrt_rev(a,b) = sqrt_rev(promote(a,b)...)
-
-
-# IEEE-1788 style
-
-function sqr_rev(c, x)   # c = x^2;  refine x
-    x1 = sqrt(c) ∩ x
-    x2 = -(sqrt(c)) ∩ x
-
-    return (c, hull(x1, x2))
-end
-
-sqr_rev(c) = sqr_rev(c, -∞..∞)
-
-"""
-∘_rev1(b, c, x) is the subset of x such that x ∘ b is defined and in c
-∘_rev2(a, c, x) is the subset of x such that a ∘ x is defined and in c
-
-If these agree (∘ is commutative) then call it ∘_rev(b, c, x)
-"""
-
-function mul_rev_new(b, c, x)   # c = b*x
-    return x ∩ (c / b)
-end
-
-function pow_rev1(b, c, x)   # c = x^b
-    return x ∩ c^(1/b)
-end
-
-function pow_rev2(a, c, x)   # c = a^x
-    return x ∩ (log(c) / lob(a))
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index bac3e49..4a9208d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -163,9 +163,3 @@ end
 
     @test C3(A, x) == IntervalBox(sqrt(A / 16))
 end
-
-@testset "power_rev for odd power" begin
-    x = -∞..∞
-    a = -8..27
-    power_rev(a, x, 3)[2] == Interval(-2.0000000000000004, 3.0000000000000004)
-end