diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 4d62b07..7fa5ebf 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -16,25 +16,13 @@ jobs:
       matrix:
 
         include:
-          - version: '1'
-            os: ubuntu-latest
-            arch: x64
-          - version: '1.0'
-            os: ubuntu-latest
-            arch: x64
-          - version: '1.1'
-            os: ubuntu-latest
-            arch: x64
-          - version: '1.2'
-            os: ubuntu-latest
-            arch: x64
-          - version: '1.3'
+          - version: '1.6'
             os: ubuntu-latest
             arch: x64
-          - version: '1.4'
-            os: ubuntu-latest
+          - version: '1.6'
+            os: windows-latest
             arch: x64
-          - version: '1.5'
+          - version: '1'
             os: ubuntu-latest
             arch: x64
           - version: '1'
diff --git a/Project.toml b/Project.toml
index 784cc5f..06c8541 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "McCormick"
 uuid = "53c679d3-6890-5091-8386-c291e8c8aaa1"
 authors = ["Matthew Wilhelm <matthew.wilhelm@uconn.edu>"]
-version = "0.11.0"
+version = "0.12"
 
 [deps]
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
@@ -11,21 +11,23 @@ IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-DiffRules = "0.0.5, 0.0.6, 0.0.7, 0.0.8, 0.0.9, 0.0.10, 0.1.0, ~1.0"
+DiffRules = "1.5"
 DocStringExtensions = "~0.8"
-ForwardDiff = "~0.5.0, ~0.6, ~0.7, ~0.8, ~0.9, ~0.10"
-IntervalArithmetic = "~0.14, ~0.15, ~0.16, ~0.17"
+ForwardDiff = "~0.10"
+IntervalArithmetic = "~0.20"
 IntervalRootFinding = "~0.5"
-NaNMath = "0.3.3"
-Reexport = "~0.2"
-StaticArrays = "~0.8, ~0.9, ~0.10, ~0.11, ~0.12"
-SpecialFunctions = "~0.8, ~0.9, ~0.10"
-julia = "~1"
+NaNMath = "0.3.5"
+NNlib = "~0.7"
+Reexport = "~0.2, ~1"
+StaticArrays = "~1.2"
+SpecialFunctions = "~1, ~2"
+julia = "~1.6, ~1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/benchmark/check_allocations.jl b/benchmark/check_allocations.jl
index 1be6f4e..9dc48e9 100644
--- a/benchmark/check_allocations.jl
+++ b/benchmark/check_allocations.jl
@@ -28,7 +28,7 @@ for T in (NS, Diff, MV)
                !iszero(allocs(@benchmark ($op)($a) samples = 1)) && println("0 in X $op Allocates")
                a = MC{5,T}(-0.5,(Interval{Float64}(-0.9,-0.1)),2)
                !iszero(allocs(@benchmark ($op)($a) samples = 1)) && println("Negative $op Allocates")
-              end
+           end
            for op in (min, max, *, -, +, /)
                x = MC{5,T}(0.4,(Interval{Float64}(0.1,0.9)),2)
                y = MC{5,T}(0.5,(Interval{Float64}(0.3,0.9)),2)
@@ -52,7 +52,6 @@ for T in (NS, Diff, MV)
            a = MC{5,T}(-0.5,Interval{Float64}(-0.9,-0.1),2); !iszero(allocs(@benchmark $a^(-2) samples = 1)) && println("a^-2 Neg Allocates")
            a = MC{5,T}(0.4,Interval{Float64}(0.1,0.9),2);    !iszero(allocs(@benchmark $a^(-3) samples = 1)) && println("a^-3 Pos Allocates")
            a = MC{5,T}(-0.5,Interval{Float64}(-0.9,-0.1),2); !iszero(allocs(@benchmark $a^(-3) samples = 1)) && println("a^-3 Neg Allocates")
-   end
 end
 
 println("Allocation Check Complete!")
diff --git a/src/McCormick.jl b/src/McCormick.jl
index b88ebae..7132af0 100644
--- a/src/McCormick.jl
+++ b/src/McCormick.jl
@@ -25,9 +25,11 @@ import Base: +, -, *, /, convert, in, isempty, one, zero, real, eps, max, min,
              cos, tan, min, max, sec, csc, cot, ^, step, sign, intersect,
              promote_rule, asinh, atanh, tanh, atan, asin, cosh, acos,
              sind, cosd, tand, asind, acosd, atand,
-             secd, cscd, cotd, asecd, acscd, acotd, isone, isnan, empty,
+             secd, cscd, cotd, asecd, acscd, acotd, isone, isnan, isfinite, empty,
              <, <=, ==, fma, cbrt, sinpi, cospi, union
 
+import NNlib: relu, selu, leakyrelu, sigmoid, swish, gelu, elu, softsign, logcosh, softplus
+
 using IntervalArithmetic
 using IntervalArithmetic: @round
 if isdefined(IntervalArithmetic, :big53)
@@ -49,12 +51,8 @@ import IntervalArithmetic: dist, mid, pow, +, -, *, /, convert, in, isempty,
                            secd, cscd, cotd, asecd, acscd, acotd, half_pi,
                            setrounding, diam, isthin, abs2
 
-if ~(VERSION < v"1.1-")
-    import SpecialFunctions: erf, erfc, erfinv, erfcinv
-    export erf, erfinv, erfc, erfcinv, erf_kernel,
-           erfinv_kernel, erfc_kernel, erfcinv_kernel
-end
-
+import SpecialFunctions: erf, erfc, erfinv, erfcinv
+export erf, erfinv, erfc, erfcinv, erf_kernel, erfinv_kernel, erfc_kernel, erfcinv_kernel
 import Base.MathConstants.golden
 
 # Export forward operators
@@ -67,7 +65,7 @@ export MC, MCNoGrad, cc, cv, Intv, lo, hi,  cc_grad, cv_grad, cnst, +, -, *, /,
        sind, cosd, tand, asind, acosd, atand, nan,
        sinhd, coshd, tanhd, asinhd, acoshd, atanhd,
        secd, cscd, cotd, asecd, acscd, acotd,
-       secdh, cschd, cothd, asechd, acschd, acothd, isone, isnan, interval_MC,
+       secdh, cschd, cothd, asechd, acschd, acothd, isone, isfinite, isnan, interval_MC,
        relu, param_relu, leaky_relu, maxsig, maxtanh, softplus, pentanh,
        sigmoid, bisigmoid, softsign, gelu, elu, selu, swish1,
        positive, negative, lower_bnd, upper_bnd, bnd, xlogx,
@@ -617,5 +615,6 @@ using Reexport
 @reexport using IntervalArithmetic
 @reexport using StaticArrays
 @reexport using SpecialFunctions
+@reexport using NNlib
 
 end
diff --git a/src/forward_operators/activation_functions.jl b/src/forward_operators/activation_functions.jl
index c24cda7..2ec6c95 100644
--- a/src/forward_operators/activation_functions.jl
+++ b/src/forward_operators/activation_functions.jl
@@ -8,16 +8,15 @@
 # src/forward_operators/concave_increasing.jl
 # Contains definitions of commonly used activation functions:
 # relu, parametric relu, leaky relu, maxsig, elu, selu, gelu, sigmoid,
-# bisigmoid, swish1, maxtanh, pentanh, softsign, softplus.
+# bisigmoid, swish, maxtanh, pentanh, softsign, softplus.
 #############################################################################
 
 # RELU DEFINITION
-"""
+#=
 relu
 
 The Rectified Linear Unit (ReLU) activation function `relu(x) = max(x, 0.0)`.
-"""
-relu(x) = max(x, 0.0)
+=#
 relu_deriv(x) = x > 0.0 ? 1.0 : 0.0
 relu_deriv2(x) =  0.0
 function relu_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
@@ -122,13 +121,11 @@ end
 maxsig(x::MC{N,T}) where {N, T<:Union{NS,MV}} = maxsig_kernel(x, maxsig(x.Intv))
 
 # DEFINE ELU
-"""
+#=
 elu
 
 The Exponential Linear Unit (ELU) activation function  `elu(x, α) = x > 0 ? x : α*(exp(x) - 1.0)`.
-"""
-@inline elu(x, α) = x > 0 ? x : α*(exp(x) - 1.0)
-@inline elu(x::Float64, α::Float64) = x > 0 ? x : α*(exp(x) - 1.0)
+=#
 function elu(x::Interval{Float64}, α::Float64)
     xL = x.lo
     xU = x.hi
@@ -220,13 +217,11 @@ function maxtanh_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,
 end
 maxtanh(x::MC{N,T}) where {N, T<:Union{NS,MV}} = maxtanh_kernel(x, maxtanh(x.Intv))
 
-"""
+#=
 softplus
 
 The `softplus` activation function  `softplus(x) = log(1.0 + exp(x))`.
-"""
-@inline softplus(x) = log(1.0 + exp(x))
-@inline softplus(x::Float64) = log(1.0 + exp(x))
+=#
 @inline softplus(x::Interval{Float64}) = log(1.0 + exp(x))
 @inline softplus_deriv(x::Float64) = 1.0/(exp(-x) + 1.0)
 @inline softplus_deriv2(x::Float64) = exp(-x)/(exp(-x) + 1.0)^2
@@ -297,18 +292,16 @@ end
     return pentanh(x), pentanh_deriv(x), p
 end
 
-"""
+#=
 sigmoid
 
 The `sigmoid` activation function `sigmoid(x) = 1.0/(1.0 + exp(-x))`.
-"""
-@inline sigmoid(x) = 1.0/(1.0 + exp(-x))
-@inline sigmoid(x::Float64) = 1.0/(1.0 + exp(-x))
+=#
 @inline sigmoid(x::Interval{Float64}) = 1.0/(1.0 + exp(-x))
-@inline sigmoid_deriv(x::Float64) = sigmoid(x)*(1.0 - sigmoid(x))
-@inline sigmoid_deriv2(x::Float64) = 2.0*exp(-2.0*x)/sigmoid(x)^3 - exp(-x)/sigmoid(x)^2
+@inline sigmoid_deriv(x::Float64) = exp(-x)/(1.0 + exp(-x))^2
+@inline sigmoid_deriv2(x::Float64) = -exp(x)*(exp(x) - 1.0)/(exp(x) + 1.0)^3
 @inline function sigmoid_env(x::Float64, y::Float64, z::Float64)
-    (x - y) - (sigmoid(x) - sigmoid(y))/sigmoid_deriv(x)
+    (x - y)*sigmoid_deriv(x) - (sigmoid(x) - sigmoid(y))
 end
 @inline function cv_sigmoid(x::Float64, xL::Float64, xU::Float64, p::Float64)
     (xL >= 0.0) && (return dline_seg(sigmoid, sigmoid_deriv, x, xL, xU)..., p)
@@ -376,13 +369,11 @@ end
     return bisigmoid(x), bisigmoid_deriv(x), p
 end
 
-"""
+#=
 softsign
 
 The `softsign` activation function `softsign(x) = x/(1.0 + abs(x))`.
-"""
-@inline softsign(x) = x/(1.0 + abs(x))
-@inline softsign(x::Float64) = x/(1.0 + abs(x))
+=#
 @inline function softsign(x::Interval{Float64})
     xLintv = Interval(x.lo)
     xUintv = Interval(x.hi)
@@ -391,16 +382,9 @@ The `softsign` activation function `softsign(x) = x/(1.0 + abs(x))`.
     return Interval(xLc.hi, xUc.hi)
 end
 @inline softsign_deriv(x::Float64) = 1.0/(1.0 + abs(x))^2
-@inline function softsign_deriv2(x::Float64)
-    if x >= 0.0
-        xp1 = 1.0 + x
-        return 2.0*(x*xp1^(-3) - xp1^(-2))
-    end
-    xm1 = 1.0 - x
-    return 2.0*(x*xm1^(-3) + xm1^(-2))
-end
+@inline softsign_deriv2(x::Float64) = -2.0*x/(abs(x)*(1.0 + abs(x))^(3))
 @inline function softsign_env(x::Float64, y::Float64, z::Float64)
-    (x - y) - (softsign(x) - softsign(y))/softsign_deriv(x)
+    (x - y)*softsign_deriv(x) - (softsign(x) - softsign(y))
 end
 @inline function cv_softsign(x::Float64, xL::Float64, xU::Float64, p::Float64)
     (xL >= 0.0) && (return dline_seg(softsign, softsign_deriv, x, xL, xU)..., p)
@@ -436,8 +420,6 @@ gelu
 
 The Gaussian Error Linear Unit `gelu` activation function `gelu(x) = x/(1.0 + abs(x))`.
 """
-@inline gelu(x) = x*(1.0 + erf(x/sqrt(2)))/2.0
-@inline gelu(x::Float64) = x*(1.0 + erf(x/sqrt(2)))/2.0
 @inline function gelu(x::Interval{Float64})
     xLintv = Interval(x.lo)
     xUintv = Interval(x.hi)
@@ -499,7 +481,7 @@ end
                 flag && (p1 = golden_section(xL, GELU_2D_ROOT1, gelu_env, xU, 0.0))
             end
             (x >= p1) && (return dline_seg(gelu, gelu_deriv, x, p1, xU)..., p1, p2)
-            return swish1(x), swish1_deriv(x), p1, p2
+            return gelu(x), gelu_deriv(x), p1, p2
         else
             return dline_seg(gelu, gelu_deriv, x, xL, xU)..., p1, p2
         end
@@ -558,13 +540,11 @@ end
 end
 
 """
-swish1
+swish
 
-The Swish-1 activation function `swish1(x) = x/(1.0 + exp(-x))`.
+The Swish-1 activation function `swish(x) = x/(1.0 + exp(-x))`.
 """
-@inline swish1(x) = x/(1.0 + exp(-x))
-@inline swish1(x::Float64) = x/(1.0 + exp(-x))
-@inline function swish1(x::Interval{Float64})
+@inline function swish(x::Interval{Float64})
     xLintv = Interval(x.lo)
     xUintv = Interval(x.hi)
     xLc = xLintv/(1.0 + exp(-xLintv))
@@ -581,40 +561,38 @@ The Swish-1 activation function `swish1(x) = x/(1.0 + exp(-x))`.
     end
     return Interval(xLcv, xUcv)
 end
-@inline function swish1_deriv(x::Float64)
-    sigmoid(x) + x*sigmoid_deriv(x)
+@inline function swish_deriv(x::Float64)
+    (x - 1)/(exp(x) + 1) - x/(exp(x) + 1)^2 + 1
 end
-@inline function swish1_deriv2(x::Float64)
-    frac1 = 2.0*exp(-2.0*x)/(exp(-x) + 1.0)^3
-    frac2 = exp(-x)/(exp(-x) + 1.0)^2
-    2.0*exp(-x)/(exp(-x) + 1.0)^2 + (frac1 - frac2)*x
+@inline function swish_deriv2(x::Float64)
+    (2 - x)/(exp(x) + 1) - 2*x/(exp(x) + 1)^3 + (3*x - 2)/(exp(x) + 1)^2
 end
-@inline function swish1_env(x::Float64, y::Float64, z::Float64)
-    swish1_deriv(x)*(y - x) + swish1(x) - swish1(y)
+@inline function swish_env(x::Float64, y::Float64, z::Float64)
+    swish_deriv(x)*(y - x) + swish(x) - swish(y)
 end
-@inline function swish1_denv(x::Float64, y::Float64, z::Float64)
-    swish1_deriv2(x)*(y - x)
+@inline function swish_denv(x::Float64, y::Float64, z::Float64)
+    swish_deriv2(x)*(y - x)
 end
 
-@inline function swish1_envm(x::Float64, y::Float64, z::Float64)
-    swish1_deriv(y)*(x - y) - (swish1(x) - swish1(y))
+@inline function swish_envm(x::Float64, y::Float64, z::Float64)
+    swish_deriv(y)*(x - y) - (swish(x) - swish(y))
 end
 
 @inline function swish_rt1(x::Float64, y::Float64, z::Float64)
-    swish1_deriv(y)*(x - y) + swish1(y) - swish1(x)
+    swish_deriv(y)*(x - y) + swish(y) - swish(x)
 end
 @inline function swish_rt1_deriv(x::Float64, y::Float64, z::Float64)
-    swish1_deriv(y) - swish1_deriv(x)
+    swish_deriv(y) - swish_deriv(x)
 end
 
-@inline function cc_swish1(x::Float64, xL::Float64, xU::Float64, p1::Float64, p2::Float64)
+@inline function cc_swish(x::Float64, xL::Float64, xU::Float64, p1::Float64, p2::Float64)
     # Single convexity regions
     if xL >= SWISH1_2D_ROOT2
-        return swish1(x), swish1_deriv(x), p1, p2
+        return swish(x), swish_deriv(x), p1, p2
     elseif xU <= SWISH1_2D_ROOT1
-        return swish1(x), swish1_deriv(x), p1, p2
+        return swish(x), swish_deriv(x), p1, p2
     elseif (SWISH1_2D_ROOT1 <= xL) && (xU <= SWISH1_2D_ROOT2)
-        return dline_seg(swish1, swish1_deriv, x, xL, xU)..., p1, p2
+        return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
     end
 
     if xL < SWISH1_2D_ROOT1
@@ -622,13 +600,13 @@ end
         flag && (p2 = golden_section(SWISH1_2D_ROOT1, 0.0, swish_rt1, xL, 0.0))
         if p2 > xU
             if p1 === Inf
-                p1, flag = secant(xL, SWISH1_2D_ROOT1, xL, SWISH1_2D_ROOT1, swish1_env, xU, 0.0)
-                flag && (p1 = golden_section(xL, SWISH1_2D_ROOT1, swish1_env, xU, 0.0))
+                p1, flag = secant(xL, SWISH1_2D_ROOT1, xL, SWISH1_2D_ROOT1, swish_env, xU, 0.0)
+                flag && (p1 = golden_section(xL, SWISH1_2D_ROOT1, swish_env, xU, 0.0))
             end
-            (x >= p1) && (return dline_seg(swish1, swish1_deriv, x, p1, xU)..., p1, p2)
-            return swish1(x), swish1_deriv(x), p1, p2
+            (x >= p1) && (return dline_seg(swish, swish_deriv, x, p1, xU)..., p1, p2)
+            return swish(x), swish_deriv(x), p1, p2
         else
-            return dline_seg(swish1, swish1_deriv, x, xL, xU)..., p1, p2
+            return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
         end
     end
 
@@ -637,31 +615,31 @@ end
         flag && (p2 = golden_section(SWISH1_MIN, SWISH1_2D_ROOT2, swish_rt1, xU, 0.0))
         if p2 < xL
             if p1 === Inf
-                p1, flag = secant(SWISH1_2D_ROOT2, xU, SWISH1_2D_ROOT2, xU, swish1_env, xL, 0.0)
-                flag && (p1 = golden_section(SWISH1_2D_ROOT2, xU, swish1_env, xL, 0.0))
+                p1, flag = secant(SWISH1_2D_ROOT2, xU, SWISH1_2D_ROOT2, xU, swish_env, xL, 0.0)
+                flag && (p1 = golden_section(SWISH1_2D_ROOT2, xU, swish_env, xL, 0.0))
             end
-            (x >= p1) && (return dline_seg(swish1, swish1_deriv, x, p1, xU)..., p1, p2)
-            return swish1(x), swish1_deriv(x), p1, p2
+            (x >= p1) && (return dline_seg(swish, swish_deriv, x, p1, xU)..., p1, p2)
+            return swish(x), swish_deriv(x), p1, p2
         else
-            return dline_seg(swish1, swish1_deriv, x, xL, xU)..., p1, p2
+            return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
         end
     end
 end
-@inline function cv_swish1(x::Float64, xL::Float64, xU::Float64, p1::Float64, p2::Float64)
+@inline function cv_swish(x::Float64, xL::Float64, xU::Float64, p1::Float64, p2::Float64)
 
     # Single convexity regions
     if xL >= SWISH1_2D_ROOT2
-        return dline_seg(swish1, swish1_deriv, x, xL, xU)..., p1, p2
+        return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
     elseif xU <= SWISH1_2D_ROOT1
-        return dline_seg(swish1, swish1_deriv, x, xL, xU)..., p1, p2
+        return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
     elseif (SWISH1_2D_ROOT1 <= xL) && (xU <= SWISH1_2D_ROOT2)
-        return swish1(x), swish1_deriv(x), p1, p2
+        return swish(x), swish_deriv(x), p1, p2
     end
 
     if xL < SWISH1_2D_ROOT1
         if p1 === Inf
-            p1, flag = newton(0.5*(SWISH1_2D_ROOT1 + SWISH1_MIN), SWISH1_2D_ROOT1, SWISH1_MIN, swish1_env, swish1_denv, xL, 0.0)
-            flag && (p1 = golden_section(SWISH1_2D_ROOT1, SWISH1_MIN, swish1_env, xL, 0.0))
+            p1, flag = newton(0.5*(SWISH1_2D_ROOT1 + SWISH1_MIN), SWISH1_2D_ROOT1, SWISH1_MIN, swish_env, swish_denv, xL, 0.0)
+            flag && (p1 = golden_section(SWISH1_2D_ROOT1, SWISH1_MIN, swish_env, xL, 0.0))
         end
     else
         p1 = -Inf
@@ -669,19 +647,19 @@ end
 
     if xU > SWISH1_2D_ROOT2
         if p2 === Inf
-            p2, flag = newton(0.5*(SWISH1_MIN + SWISH1_2D_ROOT2), SWISH1_MIN, SWISH1_2D_ROOT2, swish1_env, swish1_denv, xU, 0.0)
-            flag && (p2 = golden_section(SWISH1_MIN, SWISH1_2D_ROOT2, swish1_env, xU, 0.0))
+            p2, flag = newton(0.5*(SWISH1_MIN + SWISH1_2D_ROOT2), SWISH1_MIN, SWISH1_2D_ROOT2, swish_env, swish_denv, xU, 0.0)
+            flag && (p2 = golden_section(SWISH1_MIN, SWISH1_2D_ROOT2, swish_env, xU, 0.0))
         end
     else
         p2 = Inf
     end
 
     if x < p1
-        return dline_seg(swish1, swish1_deriv, x, xL, p1)..., p1, p2
+        return dline_seg(swish, swish_deriv, x, xL, p1)..., p1, p2
     elseif x > p2
-        return dline_seg(swish1, swish1_deriv, x, p2, xU)..., p1, p2
+        return dline_seg(swish, swish_deriv, x, p2, xU)..., p1, p2
     end
-    return swish1(x), swish1_deriv(x), p1, p2
+    return swish(x), swish_deriv(x), p1, p2
 end
 
 # define kernel and operator for sigmoid, bisigmoid, softsign, gelu
@@ -710,13 +688,13 @@ for expri in (:pentanh, :sigmoid, :bisigmoid, :softsign)
     end
 end
 
-for expri in (:swish1, :gelu)
+for expri in (:swish, :gelu)
     expri_cv = Symbol("cv_"*String(expri))
     expri_cc = Symbol("cc_"*String(expri))
     expri_kernel = Symbol(String(expri)*"_kernel")
-    if expri == swish1
+    if expri == swish
         eps_min = :(SWISH1_MIN < xL ? xL : (SWISH1_MIN > xU ? xU : SWISH1_MIN))
-        eps_max = :(swish1(xL) < swish1(xU) ? xU : xL)
+        eps_max = :(swish(xL) < swish(xU) ? xU : xL)
     else
         eps_min = :(GELU_MIN > xU ? xU : (GELU_MIN < xL ? xL : GELU_MIN))
         eps_max = :(gelu(xL) < gelu(xU) ? xU : xL)
@@ -746,7 +724,6 @@ logcosh
 
 The function `logcosh` is defined as `logcosh(x) = log(cosh(x))`.
 """
-logcosh(x::Float64) = log(cosh(x))
 logcosh_deriv(x::Float64) = tanh(x)
 logcosh_deriv2(x::Float64) = sech(x)^2
 cc_logcosh(x::Float64, xL::Float64, xU::Float64) = dline_seg(logcosh, logcosh_deriv, x, xL, xU)
diff --git a/src/forward_operators/apriori_mult.jl b/src/forward_operators/apriori_mult.jl
index 2e84356..6de8705 100644
--- a/src/forward_operators/apriori_mult.jl
+++ b/src/forward_operators/apriori_mult.jl
@@ -20,114 +20,87 @@ macro temp_grad_max(c1, x1, y1, c2, x2, y2)
         end)
 end
 
-min9(args...) = findmin(args)
-max9(args...) = findmax(args)
+min4(args...) = findmin(args)
+max4(args...) = findmax(args)
 
-@inline function mult_apriori_kernel(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64},
-                                    u1::Float64, u2::Float64, a1::Float64, a2::Float64,
-                                    v1::Float64, v2::Float64, b1::Float64, b2::Float64,
-                                    u1grad::SVector{N,Float64}, u2grad::SVector{N,Float64},
-                                    v1grad::SVector{N,Float64}, v2grad::SVector{N,Float64}) where {N, T<:Union{NS, MV}}
+const APRIORI_DELTA = 1E-6
+
+function mult_apriori_kernel(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64},
+                             u1::Float64, u2::Float64, a1::Float64, a2::Float64,
+                             u1grad::SVector{N,Float64}, u2grad::SVector{N,Float64}) where {N, T<:Union{NS, MV}}
 
     x1L = lo(x1);  x1U = hi(x1);  x2L = lo(x2);  x2U = hi(x2)
-    x1cv = x1.cv;  x1cc = x1.cc;  x2cv = x2.cv; x2cc = x2.cc
+    x1Li = Interval(lo(x1));  x1Ui = Interval(hi(x1));  x2Li = Interval(lo(x2));  x2Ui = Interval(hi(x2))
+    u1i = Interval(u1); u2i = Interval(u2)
+    a1i = Interval(a1); a2i = Interval(a2)
+    x1cv = x1.cv
+    x1cc = x1.cc
+    x2cv = x2.cv 
+    x2cc = x2.cc
+    x1i = Interval(x1.cv, x1.cc)
+    x2i = Interval(x2.cv, x2.cc)
     x1cv_grad = x1.cv_grad;  x1cc_grad = x1.cc_grad;
     x2cv_grad = x2.cv_grad;  x2cc_grad = x2.cc_grad
 
-    w = x1*x2
-    w0cv = w.cv
-    w0cc = w.cc
+    a2x1i = a2i*x1i
+    a1x2i = a1i*x2i
+    x1Ux2i = x1Ui*x2i
+    x2Ux1i = x2Ui*x1i
+    x1Lx2i = x1Li*x2i
+    x2Lx1i = x2Li*x1i
 
     # additional constraints from underestimators (convex)
-    w1cv = (x2U - a2)*u1 + (x1U - a1)*u2 + min(a2*x1cv, a2*x1cc) + min(a1*x2cv, a1*x2cc) + a1*a2 - a1*x2U - a2*x1U # GOOD
-    w2cv = (x2U - x2L)*u1 + min(x2L*x1cv, x2L*x1cc) + min(a1*x2cv, a1*x2cc) - a1*x2U #BAD
-    w3cv = (x1U - x1L)*u2 + min(a2*x1cv, a2*x1cc) + min(x1L*x2cv, x1L*x2cc) - x1U*a2
-    w4cv = (a2 - x2L)*u1 + (a1 - x1L)*u2 + min(x2L*x1cv, x2L*x1cc) + min(x1L*x2cv, x1L*x2cc) - a1*a2
-
-    # additional constraints from underestimators (concave)
-    w1cc = (x2L - a2)*u1 + (a1 - x1U)*u2 + max(a2*x1cv, a2*x1cc) + max(x1U*x2cv, x1U*x2cc) - a1*x2L
-    w2cc = (x2L - x2U)*u1 + max(a1*x2cv, a1*x2cc) + max(x2U*x1cv, x2U*x1cc) - x2L*a1
-    w3cc = (x1L - x1U)*u2 + max(a2*x1cv, a2*x1cc) + max(x1U*x2cv, x1U*x2cc) - x1L*a2
-    w4cc = (a2 - x2U)*u1 + (x1L - a1)*u2 + max(x2U*x1cv, x2U*x1cc) + max(a1*x2cv, a1*x2cc) - x1L*a2
+    w1cvi = (x2Ui - a2i)*u1i + (x1Ui - a1i)*u2i + lo(a2x1i) + lo(a1x2i) + a1i*a2i - a1i*x2Ui - a2i*x1Ui - APRIORI_DELTA
+    w2cvi = (x2Ui - x2Li)*u1i + lo(x2Lx1i) + lo(a1x2i) - a1i*x2Ui - APRIORI_DELTA
+    w3cvi = (x1Ui - x1Li)*u2i + lo(a2x1i) + lo(x1Lx2i) - x1Ui*a2i - APRIORI_DELTA
+    w4cvi = (a2i - x2Li)*u1i + (a1i - x1Li)*u2i + lo(x2Lx1i) + lo(x1Lx2i) - a1i*a2i - APRIORI_DELTA
 
-    # additional constraints from overestimators (convex)
-    # 5 fails,
-    w5cv = (x2L - b2)*v1 + (x1L - b1)*v2 + min(b2*x1cv, b2*x1cc) + min(x1L*x2cv, x1L*x2cc) + b1*b2 - b1*x2L - x1L*b2
-    w6cv = (x2L - x2U)*v1 + min(x2U*x1cv, x2U*x1cc) + min(b1*x2cv, b1*x2cc) - b1*x2L
-    w7cv = (x1L - x1U)*v2 + min(b2*x1cv, b2*x1cc) + min(x1U*x2cv, x1U*x2cc) - b2*x1L
-    w8cv = (b2 - x2U)*v1 + (b1 - x1U)*v2 + min(x2U*x1cv, x2U*x1cc) + min(x1U*x2cv, x1U*x2cc) - b1*b2
+    w1cv = lo(w1cvi)
+    w2cv = lo(w2cvi)
+    w3cv = lo(w3cvi)
+    w4cv = lo(w4cvi)
 
-    # additional constraints from overestimators (concave)
-    w5cc = (x2U - b2)*v1 + (b1 - x1L)*v2 + max(b2*x1cv, b2*x1cc) + max(x1L*x2cv, x1L*x2cc) - b1*x2U
-    w6cc = (x2U - x2L)*v1 + max(b1*x2cv, b1*x2cc) + max(x2L*x1cv, x2L*x1cc) - x1U*b1
-    w7cc = (x1U - x1L)*v2 + max(b2*x1cv, b2*x1cc) + max(x1L*x2cv, x1L*x2cc) - x1U*b2
-    w8cc = (b2 - x2L)*v1 + (x1U - b1)*v2 + max(x2L*x1cv, x2L*x1cc) + max(b1*x2cv, b1*x2cc) - x1U*b2
-
-    cv, cvind = max9(w0cv, w1cv, w2cv, w3cv, w4cv, w6cv, w7cv, w8cv)
-    #cv, cvind = max9(w0cv, w1cv, w2cv, w3cv, w4cv, w5cv, w6cv, w7cv, w8cv)
+    # additional constraints from underestimators (concave)
+    w1cci = (x2Li - a2i)*u1i + (a1i - x1Ui)*u2i + hi(a2x1i) + hi(x1Ux2i) - a1i*x2Li + APRIORI_DELTA
+    w2cci = (x2Li - x2Ui)*u1i + hi(a1x2i) + hi(x2Ux1i) - x2Li*a1i + APRIORI_DELTA
+    w3cci = (x1Li - x1Ui)*u2i + hi(a2x1i) + hi(x1Ux2i) - x1Li*a2i + APRIORI_DELTA
+    w4cci = (a2i - x2Ui)*u1i + (x1Li - a1i)*u2i + hi(x2Ux1i) + hi(a1x2i) - x1Li*a2i + APRIORI_DELTA
 
-    #cv = w8cv
-    #cvind = 9
+    w1cc = hi(w1cci)
+    w2cc = hi(w2cci)
+    w3cc = hi(w3cci)
+    w4cc = hi(w4cci)
 
+    cv, cvind = max4(w1cv, w2cv, w3cv, w4cv)
     if cvind == 1
-        cv_grad = w.cv_grad
-    elseif cvind == 2
         @temp_grad_min(a2, x1cv, x1cc, a1, x2cv, x2cc)
-        cv_grad = (x2U - a2)*u2grad + (x1U - a1)*u2grad + grad_temp
-    elseif cvind == 3
+        cv_grad = (x2U - a2)*u1grad + (x1U - a1)*u2grad + grad_temp
+    elseif cvind == 2
         @temp_grad_min(x2L, x1cv, x1cc, a1, x2cv, x2cc)
         cv_grad = (x2U - x2L)*u2grad + grad_temp
-    elseif cvind == 4
+    elseif cvind == 3
         @temp_grad_min(a2, x1cv, x1cc, x1L, x2cv, x2cc)
         cv_grad = (x1U - x1L)*u2grad + grad_temp
-    elseif cvind == 5
+    elseif cvind == 4
         @temp_grad_min(x2L, x1cv, x1cc, x1L, x2cv, x2cc)
         cv_grad = (a2 - x2L)*u1grad + (a1 - x1L)*u2grad + grad_temp
-    elseif cvind == 6
-        @temp_grad_min(b2, x1cv, x1cc, x1L, x2cv, x2cc)
-        cv_grad = (x2L - b2)*v1grad + (x1L - b1)*v2grad + grad_temp
-    elseif cvind == 7
-        @temp_grad_min(x2U, x1cv, x1cc, b1, x2cv, x2cc)
-        cv_grad = (x2L - x2U)*v1grad + grad_temp
-    elseif cvind == 8
-        @temp_grad_min(b2, x1cv, x1cc, b1, x2cv, x2cc)
-        cv_grad = (x1L - x1U)*v2grad + grad_temp
-    elseif cvind == 9
-        @temp_grad_min(x2U, x1cv, x1cc, x1U, x2cv, x2cc)
-        cv_grad = (b2 - x2U)*v1grad + (b1 - x1U)*v2grad + grad_temp
     end
 
-    # w5cc
-    cc, ccind = min9(w0cc, w1cc, w2cc, w3cc, w4cc, w6cc, w7cc, w8cc)
-
+    cc, ccind = min4(w1cc, w2cc, w3cc, w4cc)
     if ccind == 1
-        cc_grad = w.cc_grad
-    elseif ccind == 2
         @temp_grad_max(a2, x1cv, x1cc, x1U, x2cv, x2cc)
         cc_grad = (x2L - a2)*u1grad + (a1 - x1U)*u2grad + grad_temp
-    elseif ccind == 3
+    elseif ccind == 2
         @temp_grad_max(a1, x1cv, x1cc, x2U, x1cv, x1cc)
         cc_grad = (x2L - x2U)*u1grad + grad_temp
-    elseif ccind == 4
+    elseif ccind == 3
         @temp_grad_max(a2, x1cv, x1cc, x1U, x2cv, x2cc)
         cc_grad = (x1L - x1U)*u2grad + grad_temp
-    elseif ccind == 5
+    elseif ccind == 4
         @temp_grad_max(x2U, x1cv, x1cc, a1, x2cv, x2cc)
         cc_grad = (a2 - x2U)*u1grad + (x1L - a1)*u2grad + grad_temp
-    elseif ccind == 6
-        @temp_grad_max(b2, x1cv, x1cc, x1L, x2cv, x2cc)
-        cc_grad = (x2U - b2)*v1grad + (b1 - x1L)*v2grad + grad_temp
-    elseif ccind == 7
-        @temp_grad_max(b2, x2cv, x2cc, x2L, x1cv, x1cc)
-        cc_grad = (x2U - x2L)*v1grad + grad_temp
-    elseif ccind == 8
-        @temp_grad_max(b2, x1cv, x1cc, x1L, x2cv, x2cc)
-        cc_grad = (x1U - x1L)*v2grad + grad_temp
-    elseif ccind == 9
-        @temp_grad_max(x2L, x1cv, x1cc, b1, x2cv, x2cc)
-        cc_grad = (b2 - x2L)*v1grad + (x1U - b1)*v2grad + grad_temp
     end
 
-    cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
-    MC{N,T}(cv, cc, z, cv_grad, cc_grad, x1.cnst && x2.cnst)
+    out = MC{N,T}(cv, cc, z, cv_grad, cc_grad, x1.cnst && x2.cnst)
+    return out
 end
diff --git a/src/forward_operators/convex_increasing.jl b/src/forward_operators/convex_increasing.jl
index 69db889..29abdf5 100644
--- a/src/forward_operators/convex_increasing.jl
+++ b/src/forward_operators/convex_increasing.jl
@@ -21,9 +21,14 @@ for opMC in (:exp, :exp2, :exp10, :expm1)
               midcc, cc_id = mid3(x.cc, x.cv, xU)
               midcv, cv_id = mid3(x.cc, x.cv, xL)
               concave = xUc
-              (xUc > xLc) && (concave = (xLc*(xU - midcc) + xUc*(midcc - xL))/del)
+              if xUc - xLc > 3.0*eps(Float64)
+                  concave = (xLc*(xU - midcc) + xUc*(midcc - xL))/del
+                  concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*(xUc - xLc)/del
+              else
+                  concave = xUc
+                  concave_grad = mid_grad(x.cc_grad, x.cv_grad, 3)
+              end
               convex = ($opMC)(midcv)
-              concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*(xUc - xLc)/del
               convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*$dop
               convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
               return MC{N, T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
diff --git a/src/forward_operators/forward.jl b/src/forward_operators/forward.jl
index c19787a..6dbd454 100644
--- a/src/forward_operators/forward.jl
+++ b/src/forward_operators/forward.jl
@@ -17,6 +17,6 @@ include("comparison.jl")
 include("trilinear.jl")
 include("apriori_mult.jl")
 
-joinpath(@__DIR__, "no_gradient", "arithmetic.jl")
-joinpath(@__DIR__, "no_gradient", "convex_increasing.jl")
-joinpath(@__DIR__, "no_gradient", "mixed_convexity.jl")
+#joinpath(@__DIR__, "no_gradient", "arithmetic.jl")
+#joinpath(@__DIR__, "no_gradient", "convex_increasing.jl")
+#joinpath(@__DIR__, "no_gradient", "mixed_convexity.jl")
diff --git a/src/forward_operators/interval_special_functions.jl b/src/forward_operators/interval_special_functions.jl
index 70c1a9b..b204b80 100644
--- a/src/forward_operators/interval_special_functions.jl
+++ b/src/forward_operators/interval_special_functions.jl
@@ -48,9 +48,6 @@ erfinv(x, r::RoundingMode{:Up}) = _erfinv(x).hi
 erfcinv(x, r::RoundingMode{:Down}) = _erfcinv(x).hi
 erfcinv(x, r::RoundingMode{:Up}) = _erfcinv(x).lo
 
-erfinv(a::BigFloat) = mid(_erfinv(a))
-erfcinv(a::BigFloat) = mid(_erfcinv(a))
-
 function _erfinv(a::T) where T
     domain = Interval{T}(-1, 1)
     a ∉ domain && return DomainError("$a is not in [-1, 1]")
diff --git a/src/forward_operators/mixed_convexity.jl b/src/forward_operators/mixed_convexity.jl
index 70dbe1c..52ad878 100644
--- a/src/forward_operators/mixed_convexity.jl
+++ b/src/forward_operators/mixed_convexity.jl
@@ -194,7 +194,7 @@ end
     return dline_seg(NaNMath.atanh, atanh_deriv, x, p, xU)..., p
 end
 
-@inline tanh_deriv(x::Float64, y::Float64, z::Float64) = sech(x)^2
+@inline tanh_deriv(x::Float64) = sech(x)^2
 @inline tanh_env(x::Float64, y::Float64, z::Float64) = (x - y) - (tanh(x) - tanh(y))/(1.0 - tanh(x)^2)
 @inline tanh_envd(x::Float64, y::Float64, z::Float64)= 2.0*tanh(x)/(1.0 - tanh(x)^2)*(tanh(x) - tanh(y))
 @inline function cv_tanh(x::Float64, xL::Float64, xU::Float64, p::Float64)
@@ -218,7 +218,7 @@ end
     return tanh(x), sech(x)^2, p
 end
 
-@inline atan_deriv(x::Float64, y::Float64, z::Float64) = 1.0/(1.0 + x^2)
+@inline atan_deriv(x::Float64) = 1.0/(1.0 + x^2)
 @inline atan_env(x::Float64, y::Float64, z::Float64) = (x - y) - (1.0 + sqr(x))*(atan(x) - atan(y))
 @inline function cv_atan(x::Float64, xL::Float64, xU::Float64, p::Float64)
     (xL >= 0.0) && (return dline_seg(atan, atan_deriv, x, xL, xU)..., p)
@@ -241,7 +241,7 @@ end
     return atan(x), 1.0/(1.0+x^2), p
 end
 
-@inline asin_deriv(x::Float64, y::Float64, z::Float64) = 1.0/NaNMath.sqrt(1.0 - x^2)
+@inline asin_deriv(x::Float64) = 1.0/NaNMath.sqrt(1.0 - x^2)
 @inline function asin_env(x::Float64, y::Float64, z::Float64)
     return (NaNMath.asin(x) - NaNMath.asin(y))*NaNMath.sqrt(1.0-x^2) - x + y
 end
@@ -266,7 +266,7 @@ end
     return dline_seg(NaNMath.asin, asin_deriv, x, p, xU)..., p
 end
 
-@inline tan_deriv(x::Float64, y::Float64, z::Float64) = sec(x)^2
+@inline tan_deriv(x::Float64) = sec(x)^2
 @inline function tan_env(x::Float64, y::Float64, z::Float64)
     return (x - y) - (NaNMath.tan(x) - NaNMath.tan(y))/(1.0 + NaNMath.tan(x)^2)
 end
@@ -294,7 +294,7 @@ end
     return dline_seg(NaNMath.tan, tan_deriv, x, p, xU)..., p
 end
 
-@inline acos_deriv(x::Float64, y::Float64, z::Float64) = -1.0/NaNMath.sqrt(1.0-x^2)
+@inline acos_deriv(x::Float64) = -1.0/NaNMath.sqrt(1.0-x^2)
 @inline acos_env(x::Float64, y::Float64, z::Float64) = -(NaNMath.acos(x) - NaNMath.acos(y))*NaNMath.sqrt(1-x^2) - x + y
 @inline function cc_acos(x::Float64, xL::Float64, xU::Float64, p::Float64)
     (xL >= 0.0) && (return dline_seg(NaNMath.acos, acos_deriv, x, xL, xU)..., p)
@@ -317,7 +317,7 @@ end
     return dline_seg(acos, acos_deriv, x, p, xU)..., p
 end
 
-@inline asinh_deriv(x::Float64, y::Float64, z::Float64) = 1.0/sqrt(1.0 + x^2)
+@inline asinh_deriv(x::Float64) = 1.0/sqrt(1.0 + x^2)
 @inline asinh_env(x::Float64, y::Float64, z::Float64) = (x - y) - sqrt(1.0 + sqr(x))*(asinh(x) - asinh(y))
 @inline function cv_asinh(x::Float64, xL::Float64, xU::Float64, p::Float64)
     (xL >= 0.0) && (return dline_seg(asinh, asinh_deriv, x, xL, xU)..., p)
diff --git a/src/forward_operators/no_gradient/activation_functions.jl b/src/forward_operators/no_gradient/activation_functions.jl
new file mode 100644
index 0000000..76d6948
--- /dev/null
+++ b/src/forward_operators/no_gradient/activation_functions.jl
@@ -0,0 +1,185 @@
+# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
+# This code is licensed under MIT license (see LICENSE.md for full details)
+#############################################################################
+# McCormick.jl
+# A McCormick operator library in Julia
+# See https://github.com/PSORLab/McCormick.jl
+#############################################################################
+# src/forward_operators/concave_increasing.jl
+# Contains definitions of commonly used activation functions:
+# relu, parametric relu, leaky relu, maxsig, elu, selu, gelu, sigmoid,
+# bisigmoid, swish, maxtanh, pentanh, softsign, softplus.
+#############################################################################
+
+# RELU DEFINITION
+relu_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64}) = max_kernel(t, x, 0.0, z)
+function relu(t::ANYRELAX, x::MCNoGrad)
+    relu_Intv = max(x.Intv, 0.0)
+    z, cvi, cci, dcv, dcc = relu_kernel(t, x, relu_Intv)
+    return z
+end
+
+@inline leaky_relu_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64}) = param_relu_kernel(t, x, 0.01, z)
+@inline leaky_relu(t::ANYRELAX, x::MCNoGrad) = leaky_relu_kernel(t, x, param_relu(x.Intv, 0.01))
+
+function param_relu_kernel(t::ANYRELAX, x::MCNoGrad, α::Float64, z::Interval{Float64})
+    @assert α >= 0.0
+    xLc = z.lo
+    xUc = z.hi
+    xL = x.Intv.lo
+    xU = x.Intv.hi
+    midcv, cvi = mid3(x.cc, x.cv, xL)
+    midcc, cci = mid3(x.cc, x.cv, xU)
+    dcv = param_relu_deriv(midcv, α)
+    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
+    u = param_relu(midcv, α)
+    o = dcc*(midcc - xL) + xLc
+    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
+end
+function param_relu(t::ANYRELAX, x::MCNoGrad, α::Float64)
+    param_relu_kernel_Intv = param_relu(x.Intv, α)
+    z, cvi, cci, dcv, dcc = param_relu_kernel(t, x, α, param_relu_kernel_Intv)
+    return z
+end
+
+function maxtanh_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64})
+    xLc = z.lo
+    xUc = z.hi
+    xL = x.Intv.lo
+    xU = x.Intv.hi
+    midcv, cvi = mid3(x.cc, x.cv, xL)
+    midcc, cci = mid3(x.cc, x.cv, xU)
+    dcv = maxtanh_deriv(midcv)
+    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
+    u = maxtanh(midcv)
+    o = dcc*(midcc - xL) + xLc
+    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
+end
+function maxtanh(t::ANYRELAX, x::MCNoGrad) 
+    maxtanh_Intv = maxtanh(x.Intv)
+    z, cvi, cci, dcv, dcc = maxtanh_kernel(t, x, maxtanh_Intv)
+    return z
+end
+
+function elu_kernel(t::ANYRELAX, x::MCNoGrad, α::Float64, z::Interval{Float64})
+    xLc = z.lo
+    xUc = z.hi
+    xL = x.Intv.lo
+    xU = x.Intv.hi
+    midcv, cvi = mid3(x.cc, x.cv, xL)
+    midcc, cci = mid3(x.cc, x.cv, xU)
+    dcv = elu_deriv(midcv, α)
+    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
+    u = elu(midcv, α)
+    o = dcc*(midcc - xL) + xLc
+    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
+end
+function elu(t::ANYRELAX, x::MCNoGrad, α::Float64) 
+    elu_Intv = elu(x.Intv, α)
+    z, cvi, cci, dcv, dcc = elu_kernel(t, x, α, elu_Intv)
+    return z
+end
+
+function softplus_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64})
+    xLc = z.lo
+    xUc = z.hi
+    xL = x.Intv.lo
+    xU = x.Intv.hi
+    midcv, cv_id = mid3(x.cc, x.cv, xL)
+    midcc, cc_id = mid3(x.cc, x.cv, xU)
+    dcv = softplus_deriv(midcv)
+    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
+    u = softplus(midcv)
+    o = dcc*(midcc - xL) + xLc
+    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
+end
+function softplus(t::ANYRELAX, x::MCNoGrad) 
+    maxtanh_Intv = softplus(x.Intv)
+    z, cvi, cci, dcv, dcc = softplus_kernel(t, x, maxtanh_Intv)
+    return z
+end
+
+function maxsig_kernel(t::ANYRELAX, x::MCNoGrad, z::Interval{Float64})
+    xLc = z.lo
+    xUc = z.hi
+    xL = x.Intv.lo
+    xU = x.Intv.hi
+    midcv, cv_id = mid3(x.cc, x.cv, xL)
+    midcc, cc_id = mid3(x.cc, x.cv, xU)
+    dcv = maxsig_deriv(midcv)
+    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
+    u = maxsig(midcv)
+    o = dcc*(midcc - xL) + xLc
+    return MCNoGrad(u, o, z, x.cnst), cvi, cci, dcv, dcc
+end
+function maxsig(t::ANYRELAX, x::MCNoGrad)
+    maxsig_Intv = maxsig(x.Intv)
+    z, cvi, cci, dcv, dcc = maxsig_kernel(t, x, maxsig_Intv)
+    return z
+end
+
+# define kernel and operator for sigmoid, bisigmoid, softsign, gelu
+for expri in (:pentanh, :sigmoid, :bisigmoid, :softsign)
+    expri_cv = Symbol("cv_"*String(expri))
+    expri_cc = Symbol("cc_"*String(expri))
+    expri_kernel = Symbol(String(expri)*"_kernel")
+    eps_min = :xL
+    eps_max = :xU
+    @eval @inline function ($expri_kernel)(t::ANYRELAX, x::MCNoGrad, y::Interval{Float64}, cv_p::Float64, cc_p::Float64)
+        xL = x.Intv.lo
+        xU = x.Intv.hi
+        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
+        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
+        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
+        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
+        return MCNoGrad(cv, cc, y, x.cnst), cv_id, cc_id, cv_p, cc_p
+    end
+    @eval @inline function ($expri)(t::ANYRELAX, x::MCNoGrad)
+        z, tp1, tp2 = ($expri_kernel)(t, x, ($expri)(x.Intv), Inf, Inf)
+        return z
+    end
+end
+
+for expri in (:swish, :gelu)
+    expri_cv = Symbol("cv_"*String(expri))
+    expri_cc = Symbol("cc_"*String(expri))
+    expri_kernel = Symbol(String(expri)*"_kernel")
+    if expri == swish
+        eps_min = :(SWISH1_MIN < xL ? xL : (SWISH1_MIN > xU ? xU : SWISH1_MIN))
+        eps_max = :(swish(xL) < swish(xU) ? xU : xL)
+    else
+        eps_min = :(GELU_MIN > xU ? xU : (GELU_MIN < xL ? xL : GELU_MIN))
+        eps_max = :(gelu(xL) < gelu(xU) ? xU : xL)
+    end
+    @eval @inline function ($expri_kernel)(t::ANYRELAX, x::MCNoGrad, y::Interval{Float64},
+                            cv_p::Float64, cc_p::Float64, cv_p2::Float64, cc_p2::Float64)
+        xL = x.Intv.lo
+        xU = x.Intv.hi
+        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
+        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
+        cv, dcv, cv_p, cv_p2 = $(expri_cv)(midcv, xL, xU, cv_p, cv_p2)
+        cc, dcc, cc_p, cc_p2 = $(expri_cc)(midcc, xL, xU, cc_p, cc_p2)
+        return MCNoGrad(cv, cc, y, x.cnst), cv_id, cc_id, cv_p, cc_p, cv_p, cc_p, cv_p2, cc_p2
+    end
+    @eval @inline function ($expri)(t::ANYRELAX, x::MCNoGrad)
+        z, cv_id, cc_id, dcv, dcc, tp1, tp2, tp3, tp4 = ($expri_kernel)(t, x, ($expri)(x.Intv), Inf, Inf, Inf, Inf)
+        return z
+    end
+end
+
+@inline function logcosh_kernel(t::ANYRELAX, x::MCNoGrad, y::Interval{Float64})
+    xL = x.Intv.lo
+    xU = x.Intv.hi
+    eps_min = in(0.0, x) ? 0.0 : (xU <= 0.0 ? xU : xL)
+    eps_max = (abs(xL) < abs(xU)) ? xU : xL
+    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
+    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
+    cv, dcv = cv_logcosh(midcv, xL, xU)
+    cc, dcc = cc_logcosh(midcc, xL, xU)
+    return MCNoGrad(cv, cc, y, x.cnst), cv_id, cc_id, dcv, dcc
+end
+function logcosh(t::ANYRELAX, x::MCNoGrad)
+    logcosh_Intv = logcosh(x.Intv)
+	z, cv_id, cc_id, dcv, dcc = logcosh_kernel(t, x, logcosh_Intv)
+    return z
+end
\ No newline at end of file
diff --git a/src/forward_operators/no_gradient/comparison.jl b/src/forward_operators/no_gradient/comparison.jl
new file mode 100644
index 0000000..169168a
--- /dev/null
+++ b/src/forward_operators/no_gradient/comparison.jl
@@ -0,0 +1,32 @@
+# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
+# This code is licensed under MIT license (see LICENSE.md for full details)
+#############################################################################
+# McCormick.jl
+# A McCormick operator library in Julia
+# See https://github.com/PSORLab/McCormick.jl
+#############################################################################
+# src/forward_operators/comparison.jl
+# Defines :<, :<=, :!=, :==.
+#############################################################################
+
+for f in (:<, :<=)
+    @eval begin
+        function ($f)(::ANYRELAX, x::MCNoGrad, y::MCNoGrad)
+            ($f)(x.cv,y.cv) && ($f)(x.cc,y.cc) && ($f)(x.Intv, y.Intv)
+        end
+        function ($f)(::ANYRELAX, x::MCNoGrad, c::Float64)
+            ($f)(x.cv, c) && ($f)(x.cc, c) && ($f)(x.Intv, c)
+        end
+        function ($f)(::ANYRELAX, c::Float64, y::MCNoGrad)
+            ($f)(c, y.cv) && ($f)(c, y.cc) && ($f)(c, y.Intv)
+        end
+    end
+end
+
+function ==(::ANYRELAX, x::MCNoGrad, y::MCNoGrad)
+    (x.cv == y.cv) && (x.cc == y.cc) && (x.Intv == y.Intv)
+end
+function ==(::ANYRELAX, x::MCNoGrad, c::Float64)
+    (x.cv == c) && (x.cc == c) && (x.Intv == c)
+end
+==(::ANYRELAX, c::Float64, y::MCNoGrad) = (y == c)
diff --git a/src/forward_operators/no_gradient/division.jl b/src/forward_operators/no_gradient/division.jl
new file mode 100644
index 0000000..876a1f0
--- /dev/null
+++ b/src/forward_operators/no_gradient/division.jl
@@ -0,0 +1,53 @@
+# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
+# This code is licensed under MIT license (see LICENSE.md for full details)
+#############################################################################
+# McCormick.jl
+# A McCormick operator library in Julia
+# See https://github.com/PSORLab/McCormick.jl
+#############################################################################
+# src/forward_operators/division.jl
+# Contains definitions of division.
+#############################################################################
+
+#=
+@inline function div_MV(::Diff, x::MCNoGrad, y::MCNoGrad, z::Interval{Float64})
+    if (0.0 < y.Intv.lo)
+        cv, cv_grad = div_diffcv(x, y)
+        cc, cc_grad = div_diffcv(-x, y)
+        return MC{N,Diff}(cv, -cc, z, cv_grad, -cc_grad, x.cnst && y.cnst)
+    end
+    cv, cv_grad = div_diffcv(-x, -y)
+    cc, cc_grad = div_diffcv(-x, y)
+    return MC{N,Diff}(cv, -cc, z, cv_grad, -cc_grad, x.cnst && y.cnst)
+end
+
+@inline function div_kernel(t::Union{NS,MV}, x::MCNoGrad, y::MCNoGrad, z::Interval{Float64})
+    if (x === y)
+        zMC = one(MCNoGrad)
+    else
+        zMC = mult_kernel(t, x, inv(y), z)
+    end
+    return zMC
+end
+
+@inline function div_kernel(t::ANYRELAX, x::MCNoGrad, y::MCNoGrad, z::Interval{Float64})
+    degen1 = (x.Intv.hi - x.Intv.lo == 0.0)
+    degen2 = (y.Intv.hi - y.Intv.lo == 0.0)
+    if x === y
+        zMC = one(MCNoGrad)
+    elseif  !degen1 || !degen2
+        zMC = div_MV(t, x, y, z)
+    else
+        zMC = mult_kernel(t, x, inv(y), z)
+    end
+    return zMC
+end
+=#
+
+@inline function /(t::ANYRELAX, x::MCNoGrad, y::MCNoGrad)
+    if 0.0 ∉ y
+        z, _ =  div_kernel(t, x, y, x.Intv/y.Intv)
+        return z
+    end
+    return nan(MCNoGrad)
+end
diff --git a/src/forward_operators/no_gradient/mixed_convexity.jl b/src/forward_operators/no_gradient/mixed_convexity.jl
index a062e1a..2483791 100644
--- a/src/forward_operators/no_gradient/mixed_convexity.jl
+++ b/src/forward_operators/no_gradient/mixed_convexity.jl
@@ -29,3 +29,62 @@ function sin(::ANYRELAX, x::MCNoGrad)
     y, _, _, _, _ = sin_kernel(x, sin(x.Intv), Inf, Inf, Inf, Inf)
     return y
 end
+
+
+# basic method overloading operator (sinh, tanh, atanh, asinh), convexoconcave or concavoconvex
+eps_min_dict = Dict{Symbol,Symbol}(:sinh => :xL, :tanh => :xL, :asinh => :xL,
+                                 :atanh => :xL, :tan => :xL, :acos => :xU,
+                                 :asin => :xL, :atan => :xL, :erf => :xL,
+                                 :cbrt => :xL, :erfinv => :xL, :erfc => :xU)
+eps_max_dict = Dict{Symbol,Symbol}(:sinh => :xU, :tanh => :xU, :asinh => :xU,
+                                 :atanh => :xU, :tan => :xU, :acos => :xL,
+                                 :asin => :xU, :atan => :xU, :erf => :xU,
+                                 :cbrt => :xU, :erfinv => :xU, :erfc => :xL)
+
+for expri in (:sinh, :tanh, :asinh, :atanh, :tan, :acos, :asin, :atan,
+              (:(SpecialFunctions.erf), :erf), :cbrt,
+              (:(SpecialFunctions.erfinv), :erfinv),
+              (:(SpecialFunctions.erfc), :erfc))
+    if expri isa Symbol
+        expri_name = expri
+        expri_sym = expri
+    else
+        expri_name = expri[1]
+        expri_sym = expri[2]
+    end
+    expri_kernel = Symbol(String(expri_sym)*"_kernel")
+    expri_cv = Symbol("cv_"*String(expri_sym))
+    expri_cc = Symbol("cc_"*String(expri_sym))
+    eps_min = eps_min_dict[expri_sym]
+    eps_max = eps_max_dict[expri_sym]
+    @eval @inline function ($expri_kernel)(t::Union{NS,MV}, x::MCNoGrad, y::Interval{Float64}, cv_p::Float64, cc_p::Float64)
+        xL = x.Intv.lo
+        xU = x.Intv.hi
+        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
+        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
+        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
+        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
+        return MCNoGrad(cv, cc, y, x.cnst), dcv, dcc, cv_p, cc_p
+    end
+    @eval @inline function ($expri_kernel)(t::Diff, x::MCNoGrad, y::Interval{Float64}, cv_p::Float64, cc_p::Float64)
+        xL = x.Intv.lo
+        xU = x.Intv.hi
+        midcv, cv_id = mid3(x.cv, x.cc, $eps_min)
+        midcc, cc_id = mid3(x.cv, x.cc, $eps_max)
+        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
+        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
+        gcv1, gdcv1, cv_p = $(expri_cv)(x.cv, xL, xU, cv_p)
+        gcc1, gdcc1, cc_p = $(expri_cc)(x.cv, xL, xU, cc_p)
+        gcv2, gdcv2, cv_p = $(expri_cv)(x.cc, xL, xU, cv_p)
+        gcc2, gdcc2, cc_p = $(expri_cc)(x.cc, xL, xU, cc_p)
+        c_cv1 = max(0.0, gdcv1)
+        c_cv2 = min(0.0, gdcv2)
+        c_cc1 = min(0.0, gdcc1)
+        c_cc2 = max(0.0, gdcc2)
+        return MCNoGrad(cv, cc, y, cv_grad, cc_grad, x.cnst), c_cv1, c_cv2, c_cc1, c_cc2
+    end
+    @eval @inline function ($expri_name)(t::ANYRELAX, x::MCNoGrad)
+        z, tp1, tp2 = ($expri_kernel)(t, x, ($expri_name)(x.Intv), Inf, Inf)
+        return z
+    end
+end
\ No newline at end of file
diff --git a/src/forward_operators/no_gradient/multiplication.jl b/src/forward_operators/no_gradient/multiplication.jl
new file mode 100644
index 0000000..27a3ebc
--- /dev/null
+++ b/src/forward_operators/no_gradient/multiplication.jl
@@ -0,0 +1,252 @@
+# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
+# This code is licensed under MIT license (see LICENSE.md for full details)
+#############################################################################
+# McCormick.jl
+# A McCormick operator library in Julia
+# See https://github.com/PSORLab/McCormick.jl
+#############################################################################
+# src/forward_operators/multiplication.jl
+# Defines multiplication operator.
+#############################################################################
+
+# Nonsmooth multiplication kernel definition
+@inline function mul1_u1pos_u2pos(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}, cnst::Bool)
+    xLc = z.lo
+    xUc = z.hi
+  	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
+  	cv2 = x2.Intv.lo*x1.cv + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
+  	if cv1 > cv2
+    	cv = cv1
+    	cv_grad = x2.Intv.hi*x1.cv_grad
+  	else
+    	cv = cv2
+    	cv_grad = x2.Intv.lo*x1.cv_grad
+  	end
+
+  	cc1 = x2.Intv.lo*x1.cc + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
+  	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
+  	if cc1 < cc2
+    	cc = cc1
+    	cc_grad = x2.Intv.lo*x1.cc_grad
+  	else
+    	cc = cc2
+    	cc_grad = x2.Intv.hi*x1.cc_grad
+  	end
+  	return MCNoGrad(cv, cc, z, cnst)
+end
+@inline function mul1_u1pos_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}, cnst::Bool)
+    xLc = z.lo
+    xUc = z.hi
+    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
+    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
+    if cv1 > cv2
+        cv = cv1
+        cv_grad = x2.Intv.hi*x1.cv_grad
+    else
+        cv = cv2
+        cv_grad = x2.Intv.lo*x1.cc_grad
+    end
+
+    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
+    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
+    if cc1 < cc2
+        cc = cc1
+        cc_grad = x2.Intv.lo*x1.cv_grad
+    else
+        cc = cc2
+        cc_grad = x2.Intv.hi*x1.cc_grad
+    end
+    return MCNoGrad(cv, cc, z, cnst)
+end
+@inline function mul1_u1mix_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}, cnst::Bool)
+    xLc = z.lo
+    xUc = z.hi
+    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
+    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.lo
+    cv = cv1
+    if cv1 > cv2
+        cv = cv1
+        cv_grad = x2.Intv.hi*x1.cv_grad
+    else
+        cv = cv2
+        cv_grad = x2.Intv.lo*x1.cc_grad
+    end
+    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
+    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.hi
+    cc = cc1
+    if cc1 < cc2
+        cc = cc1
+        cc_grad = x2.Intv.lo*x1.cv_grad
+    else
+        cc = cc2
+        cc_grad = x2.Intv.hi*x1.cc_grad
+    end
+    return MCNoGrad(cv, cc, z, cnst)
+end
+@inline function mul2_u1pos_u2pos(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}) 
+	xLc = z.lo
+	xUc = z.hi
+	cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
+	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
+	cv2 = x2.Intv.lo*x1.cv + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
+	if cv1 > cv2
+		cv = cv1
+		cv_grad = x2.Intv.hi*x1.cv_grad + x1.Intv.hi*x2.cv_grad
+	else
+		cv = cv2
+		cv_grad = x2.Intv.lo*x1.cv_grad + x1.Intv.lo*x2.cv_grad
+	end
+	cc1 = x2.Intv.lo*x1.cc + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
+	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
+	if cc1 < cc2
+		cc = cc1
+		cc_grad = x2.Intv.lo*x1.cc_grad + x1.Intv.hi*x2.cc_grad
+	else
+		cc = cc2
+		cc_grad = x2.Intv.hi*x1.cc_grad + x1.Intv.lo*x2.cc_grad
+	end
+	return MCNoGrad(cv, cc, z, cnst)
+end
+@inline function mul2_u1pos_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64}, cnst::Bool)
+    xLc = z.lo
+    xUc = z.hi
+    cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
+    cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
+    if cv1 > cv2
+        cv = cv1
+		cvt2 = x1.Intv.hi
+		cv2_is_cvg = 1
+    else
+        cv = cv2
+		cvt2 = x1.Intv.lo
+		cv2_is_cvg = 0
+    end
+	cvt1 = 0.0
+	cv1_is_cvg = 2
+
+    cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
+    cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
+    if cc1 < cc2
+        cc = cc1
+		cct2 = x1.Intv.hi
+    else
+        cc = cc2
+		cct2 = x1.Intv.lo
+    end
+	cct1 = 0.0
+	cc1_is_ccg = 2
+	cc2_is_ccg = 1
+
+	return MCNoGrad(cv, cc, z, cnst), cvt1, cvt2, cct1, cct2, cv1_is_cvg, cv2_is_cvg, cc1_is_ccg, cc2_is_ccg
+end
+@inline function mul2_u1mix_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64})
+	xLc = z.lo
+	xUc = z.hi
+  	cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
+	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
+	cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.lo
+	if cv1 > cv2
+		cv = cv1
+		cvt1 = x2.Intv.hi
+		cvt2 = x1.Intv.hi
+		cv1_is_cvg = 1
+		cv2_is_cvg = 1
+	else
+		cv = cv2
+		cvt1 = x2.Intv.lo
+		cvt2 = x1.Intv.lo
+		cv1_is_cvg = 0
+		cv2_is_cvg = 0
+	end
+
+	cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
+	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.hi
+	if cc1 < cc2
+		cc = cc1
+		cct1 = x2.Intv.lo
+		cct2 = x1.Intv.hi
+		cc1_is_ccg = 0
+		cc2_is_ccg = 1
+	else
+		cc = cc2
+		cct1 = x2.Intv.hi
+		cct2 = x1.Intv.lo
+		cc1_is_ccg = 1
+		cc2_is_ccg = 0
+	end
+
+	return MCNoGrad(cv, cc, z, cnst), cvt1, cvt2, cct1, cct2, cv1_is_cvg, cv2_is_cvg, cc1_is_ccg, cc2_is_ccg
+end
+@inline function mul3_u1pos_u2mix(t::NS, x1::MCNoGrad, x2::MCNoGrad, z::Interval{Float64})
+	xLc = z.lo
+	xUc = z.hi
+    cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
+	cv1 = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
+	cv2 = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
+	if cv1 > cv2
+		cv = cv1
+		cvt1 = x2.Intv.hi
+		cvt2 = x1.Intv.hi
+		cv1_is_cvg = 1
+	else
+		cv = cv2
+		cvt1 = x2.Intv.lo
+		cvt2 = x1.Intv.lo
+		cv1_is_cvg = 0
+	end
+	cv2_is_cvg = 1
+
+	cc1 = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
+	cc2 = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
+	if cc1 < cc2
+		cc = cc1
+		cct1 = x2.Intv.lo
+		cct2 = x1.Intv.hi
+		cc1_is_ccg = 0
+	else
+		cc = cc2
+		cct1 = x2.Intv.hi
+		cct2 = x1.Intv.lo
+		cc1_is_ccg = 1
+	end
+	cc2_is_ccg = 1
+	return MCNoGrad(cv, cc, z, cnst), cvt1, cvt2, cct1, cct2, cv1_is_cvg, cv2_is_cvg, cc1_is_ccg, cc2_is_ccg
+end
+
+@inline function multiply_STD_NS(t::NS, x1::MCNoGrad, x2::MCNoGrad, y::Interval{Float64})
+	if x2.Intv.lo >= 0.0
+    	x2.cnst && (return mul1_u1pos_u2pos(t, x1, x2, y, x1.cnst))
+    	x1.cnst && (return mul1_u1pos_u2pos(t, x2, x1, y, x2.cnst))
+  		return mul2_u1pos_u2pos(t, x1, x2, y)
+	elseif x2.Intv.hi <= 0.0
+		return -mult_kernel(t, x1, -x2, -y)
+	else
+    	x2.cnst && (return mul1_u1pos_u2mix(t, x1, x2, y, x1.cnst))
+    	x1.cnst && (return mul2_u1pos_u2mix(t, x1, x2, y, x2.cnst))
+	end
+	mul3_u1pos_u2mix(t, x1, x2, y)
+end
+
+@inline function mult_kernel(t::NS, x1::MCNoGrad, x2::MCNoGrad, y::Interval{Float64})
+	isone(x1) && (return x2)
+	isone(x2) && (return x1)
+	if x1.Intv.lo >= 0.0
+		return multiply_STD_NS(t, x1, x2, y)
+	elseif x1.Intv.hi <= 0.0
+		(x2.Intv.lo >= 0.0) && (return -mult_kernel(t, -x1, x2, -y))
+	    (x2.Intv.hi <= 0.0) && (return mult_kernel(t, -x1, -x2, y))
+		return -mult_kernel(t, x2, -x1, -y)
+	elseif x2.Intv.lo >= 0.0
+		return mult_kernel(t, x2, x1, y)
+	elseif x2.Intv.hi <= 0.0
+		return -mult_kernel(t, -x2, x1, -y)
+	else
+    	x2.cnst && (return mul1_u1mix_u2mix(t, x1, x2, y, x1.cnst))
+    	x1.cnst && (return mul1_u1mix_u2mix(t, x2, x1, y, x2.cnst))
+	end
+	return mul2_u1mix_u2mix(t, x1, x2, y)
+end
+
+@inline function *(t::NS, x1::MCNoGrad, x2::MCNoGrad)
+	mult_kernel(t, x1, x2, x1.Intv*x2.Intv)
+end
\ No newline at end of file
diff --git a/src/forward_operators/no_gradient/other.jl b/src/forward_operators/no_gradient/other.jl
new file mode 100644
index 0000000..c02d849
--- /dev/null
+++ b/src/forward_operators/no_gradient/other.jl
@@ -0,0 +1,315 @@
+# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
+# This code is licensed under MIT license (see LICENSE.md for full details)
+#############################################################################
+# McCormick.jl
+# A McCormick operator library in Julia
+# See https://github.com/PSORLab/McCormick.jl
+#############################################################################
+# src/forward_operators/other.jl
+# Defines isempty, empty, isnan, step, sign, abs, intersect, in.
+#############################################################################
+
+
+empty(x::MCNoGrad) = MCNoGrad(Inf, -Inf, Interval{Float64}(Inf,-Inf), false)
+interval_MC(x::MCNoGrad) = MCNoGrad(x.Intv)
+
+isnan(x::MCNoGrad) = isnan(x.cc) || isnan(x.cv)
+isfinite(x::MCNoGrad) = isfinite(x.cc) && isfinite(x.cv) && isfinite(x.Intv)
+
+@inline function step_kernel(t::Union{NS,MV}, x::MCNoGrad, z::Interval{Float64})
+	xL = x.Intv.lo
+	xU = x.Intv.hi
+	midcc, cc_id = mid3(x.cc, x.cv, xU)
+	midcv, cv_id = mid3(x.cc, x.cv, xL)
+	cc, dcc = cc_step_NS(midcc, xL, xU)
+	cv, dcv = cv_step_NS(midcv, xL, xU)
+	return MCNoGrad(cv, cc, z, x.cnst), cv_id, cc_id, dcv, dcc
+end
+@inline step(t::Union{NS,MV}, x::MCNoGrad) = step_kernel(t, x, step(x.Intv))
+
+@inline function abs_kernel(t::Union{NS,MV}, x::MC{N, T}, z::Interval{Float64})
+	xL = x.Intv.lo
+	xU = x.Intv.hi
+	eps_min = mid3v(xL, x.Intv.hi, 0.0)
+	eps_max = (abs(xU) >= abs(xL)) ? xU : xL
+	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
+	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
+	cc, dcc = cc_abs(midcc, xL, xU)
+	cv, dcv = cv_abs_NS(midcv, xL, xU)
+	return MCNoGrad(cv, cc, z, x.cnst), dcv, dcc
+end
+@inline abs(t::Union{NS,MV}, x::MCNoGrad) = abs_kernel(t, x, abs(x.Intv))
+
+@inline abs2_kernel(t::Union{NS,MV}, x::MCNoGrad, z::Interval{Float64}) = sqr_kernel(t, x, z)
+@inline abs2(t::Union{NS,MV}, x::MCNoGrad) = abs2_kernel(t, x, pow(x.Intv,2))
+
+@inline function xlogx_kernel(t::Union{NS,MV}, x::MCNoGrad, y::Interval{Float64})
+    xL = x.Intv.lo
+    xU = x.Intv.hi
+    eps_max = xlogx(xU) > xlogx(xL) ?  xU : xL
+    eps_min = in(1.0/exp(1.0), x) ? 1.0/exp(1.0) : (xlogx(xU) > xlogx(xL) ?  xL : xU)
+    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
+    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
+    cc, dcc = cc_xlogx(midcc, xL, xU)
+    cv, dcv = cv_xlogx(midcv, xL, xU)
+    return MCNoGrad(cv, cc, y, x.cnst), cv_id, cc_id, dcv, dcc
+end
+xlogx(t::Union{NS,MV}, x::MCNoGrad) = xlogx_kernel(t, x, xlogx(x.Intv))
+
+function arh_kernel(t::Union{NS,MV}, x::MCNoGrad, k::Float64, z::Interval{Float64}, cv_p::Float64, cc_p::Float64)
+	(k == 0.0) && return one(MC{N,T})
+	in(0.0, x) && throw(DomainError(0.0))
+	xL = x.Intv.lo
+    xU = x.Intv.hi
+	eps_min = xL
+	eps_max = xU
+    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
+    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
+    cv, dcv, cv_p = cv_arh(midcv, xL, xU, k, cv_p)
+    cc, dcc, cc_p = cc_arh(midcc, xL, xU, k, cc_p)
+    return MCNoGrad(cv, cc, z, x.cnst), cv_id, cc_id, cv_p, cc_p
+end
+function arh_kernel(x::Float64, k::MC{N,T}, z::Interval{Float64}, cv_p::Float64, cc_p::Float64) where {N,T<:Union{NS,MV}}
+	zMC, cv_id, cc_id = exp(-k/x)
+	zMC, cv_id, cc_id, 0.0, 0.0
+end
+function arh(x::MC{N,T}, k::Float64) where {N, T <: RelaxTag}
+	yMC, _, _, _, _ = arh_kernel(x, k, exp(-k/x.Intv), Inf, Inf)
+	return yMC
+end
+
+function xexpax_kernel(t::Union{NS,MV}, x::Float64, a::MCNoGrad, z::Interval{Float64}, cv_p::Float64, cc_p::Float64)
+	# TODO...
+	v1 = *(t, a, x)
+	v2 = exp(t, v1)
+	v3 = x*v2
+	v3, 0, 0, 0.0, 0.0
+end
+function xexpax_kernel(t::Union{NS,MV}, x::MCNoGrad, a::Float64, z::Interval{Float64}, cv_p::Float64, cc_p::Float64)
+	if a == 0.0
+		return one(MCNoGrad), 0, 0, 0.0, 0.0
+	end
+	xL = x.Intv.lo
+    xU = x.Intv.hi
+	zpnt = -1.0/a
+	fxL = xexpax(xL, a)
+	fxU = xexpax(xU, a)
+	if a > 0.0
+		eps_min = (xL <= zpnt <= xU) ? zpnt : ((fxL <= fxU) ? xL : xU)
+		eps_max = (fxL >= fxU) ? xL : xU
+	else
+		eps_min = (fxL <= fxU) ? xL : xU
+		eps_max = (xL <= zpnt <= xU) ? zpnt : ((fxL >= fxU) ? xL : xU)
+	end
+    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
+    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
+    cv, dcv, cv_p = cv_xexpax(midcv, xL, xU, a, cv_p)
+    cc, dcc, cc_p = cc_xexpax(midcc, xL, xU, a, cc_p)
+    return MCNoGrad(cv, cc, z, x.cnst), cv_id, cc_id, cv_p, cc_p
+end
+xexpax(t::Union{NS,MV}, x::Float64, a::MCNoGrad) = x*exp(a*x)
+function xexpax(t::Union{NS,MV}, x::MCNoGrad, a::Float64)
+	intv_xexpax = xexpax(x.Intv, a)
+	yMC, _, _, _, _ = xexpax_kernel(t, x, a, intv_xexpax, Inf, Inf)
+	return yMC
+end
+
+#=
+########### Defines sign
+@inline function correct_intersect(x::MC{N,T}, cv::Float64, cc::Float64, Intv::Interval{Float64}, cv_grad::SVector{N,Float64},
+	                               cc_grad::SVector{N,Float64}, cnst::Bool) where {N, T <: RelaxTag}
+
+    if cv - cc < MC_INTERSECT_TOL
+		if diam(Intv) > 3.0*MC_INTERSECT_TOL
+			cv = max(cv - MC_INTERSECT_TOL, Intv.lo)
+			cc = min(cc + MC_INTERSECT_TOL, Intv.hi)
+			if cv === Intv.lo
+				cv_grad  = zero(SVector{N,Float64})
+			end
+			if cc === Intv.hi
+				cc_grad = zero(SVector{N,Float64})
+			end
+			return MC{N,T}(cv, cc, Intv, cv_grad, cc_grad, cnst)
+		else
+			MC_INTERSECT_NOOP_FALLBACK && (return x)
+		end
+	end
+	return MC{N,T}(NaN, NaN, Intv, cv_grad, cc_grad, false)
+end
+
+@inline function intersect(x::MC{N,T}, y::MC{N,T}) where {N, T<:Union{NS, MV}}
+
+     Intv = x.Intv ∩ y.Intv
+	 isempty(Intv) && (return empty(x))
+
+	 cnst1 = (x.cc < y.cc) ? x.cnst : y.cnst
+	 cnst2 = (x.cv > y.cv) ? x.cnst : y.cnst
+
+     cc = (x.cc < y.cc) ? x.cc : y.cc
+     cc_grad = (x.cc < y.cc) ? x.cc_grad : y.cc_grad
+
+     cv = (x.cv > y.cv) ? x.cv : y.cv
+     cv_grad = (x.cv > y.cv) ? x.cv_grad : y.cv_grad
+
+	 if cv <= cc
+		 return MC{N,T}(cv, cc, Intv, cv_grad, cc_grad, cnst1 && cnst2)
+	 end
+
+	 return correct_intersect(x, cv, cc, Intv, cv_grad, cc_grad, cnst1 && cnst2)
+end
+
+@inline function intersect(x::MC{N, Diff}, y::MC{N, Diff}) where N
+
+	Intv = intersect(x.Intv, y.Intv)
+	isempty(Intv) && (return empty(x))
+
+    max_MC = x - max(x - y, 0.0)
+    min_MC = y - max(y - x, 0.0)
+
+	if max_MC.cv <= min_MC.cc
+    	return MC{N, Diff}(max_MC.cv, min_MC.cc, Intv, max_MC.cv_grad, min_MC.cc_grad, x.cnst && y.cnst)
+	end
+
+	return correct_intersect(x, max_MC.cv, min_MC.cc, Intv, max_MC.cv_grad, min_MC.cc_grad, x.cnst && y.cnst)
+end
+
+@inline function intersect(x::MC{N, T}, y::Interval{Float64}) where {N, T<:Union{NS,MV}}
+	cnst1 = x.cnst
+	cnst2 = x.cnst
+	if x.cv >= y.lo
+  		cv = x.cv
+  		cv_grad = x.cv_grad
+	else
+		cnst1 = true
+  		cv = y.lo
+  		cv_grad = zero(SVector{N,Float64})
+	end
+	if x.cc <= y.hi
+  		cc = x.cc
+  		cc_grad = x.cc_grad
+	else
+		cnst2 = true
+  		cc = y.hi
+  		cc_grad = zero(SVector{N,Float64})
+	end
+	if cv <= cc
+		return MC{N, T}(cv, cc, intersect(x.Intv, y), cv_grad, cc_grad, cnst1 && cnst2)
+	end
+	return correct_intersect(x, cv, cc, intersect(x.Intv, y), cv_grad, cc_grad, cnst1 && cnst2)
+end
+@inline function intersect(x::MC{N, Diff}, y::Interval{Float64}) where N
+     max_MC = x - max(x - y, 0.0)
+     min_MC = y - max(y - x, 0.0)
+	 if max_MC.cv <= min_MC.cc
+		 MC{N, Diff}(max_MC.cv, min_MC.cc, intersect(x.Intv, y), max_MC.cv_grad, min_MC.cc_grad, x.cnst)
+	 end
+     return correct_intersect(x, max_MC.cv, min_MC.cc, intersect(x.Intv, y), max_MC.cv_grad, min_MC.cc_grad, x.cnst && y.cnst)
+end
+
+@inline function intersect(c::Float64, x::MC{N,T}) where {N, T<:RelaxTag}
+	isempty(x) && (return empty(x))
+	isnan(x) && (return nan(x))
+
+	intv = intersect(x.Intv, Interval{Float64}(c))
+	isempty(intv) && (return empty(x))
+
+	cv = max(c, x.cv)
+	cc = min(c, x.cc)
+
+	cv_grad = (cv == c) ? zero(SVector{N,Float64}) : x.cv_grad
+	cc_grad = (cc == c) ? zero(SVector{N,Float64}) : x.cc_grad
+	if cv <= c <= cc
+		return MC{N, T}(c, c, intv, cv_grad, cc_grad, true)
+	end
+	correct_intersect(x, cv, cc, intv, cv_grad, cc_grad, cnst1 && cnst2)
+end
+@inline intersect(x::MC{N,T}, c::Float64) where {N, T<:RelaxTag} = intersect(c, x)
+@inline intersect(c::Interval{Float64}, x::MC{N,T}) where {N, T<:RelaxTag} = intersect(c, x)
+
+@inline in(a::Int, x::MC) = in(a, x.Intv)
+@inline in(a::T, x::MC) where T<:AbstractFloat = in(a, x.Intv)
+@inline isempty(x::MC) = isempty(x.Intv)
+
+"""
+expax
+
+The `expax` function is defined as `expax(x, a) = x*exp(a*x)`.
+
+Form defined in Najman, Jaromił, Dominik Bongartz, and Alexander Mitsos.
+"Relaxations of thermodynamic property and costing models in process engineering."
+Computers & Chemical Engineering 130 (2019): 106571.
+"""
+
+
+"""
+xabsx
+
+The function `xabsx` is defined as `xabsx(x) = x*abs(x)`.
+"""
+xabsx(x::Float64) = x*abs(x)
+xabsx_deriv(x::Float64) = 2*abs(x)
+xabsx_deriv2(x::Float64) = 2*abs(x)/x
+@inline function cc_xabsx(x::Float64, xL::Float64, xU::Float64)
+	(xU <= 0.0) && (return xabsx(x), xabsx_deriv(x))
+	(0.0 <= xL) && (return dline_seg(xabsx, xabsx_deriv, x, xL, xU))
+
+	root_f(x, xL, xU) = (xabsx(xU) - xabsx(x)) - 2*(xU - x)*abs(x)
+	root_df(x, xL, xU) = xabsx_deriv(x) - xU*xabsx_deriv2(x)
+	inflection, flag = newton(xL, xL, 0.0, root_f, root_df, xL, xU)
+	flag && (inflection = golden_section(xL, xU, root_f, xL, xU))
+	if x <= inflection
+		return xabsx(x), xabsx_deriv(x)
+	else
+		return dline_seg(xabsx, xabsx_deriv, x, inflection, xU)
+	end
+end
+@inline function cv_xabsx(x::Float64, xL::Float64, xU::Float64)
+	(xU <= 0.0) && (return dline_seg(xabsx, xabsx_deriv, x, xL, xU))
+	(0.0 <= xL) && (return xabsx(x), xabsx_deriv(x))
+
+	root_f(x, xL, xU) = (xabsx(x) - xabsx(xL)) - 2*(x - xL)*abs(x)
+	root_df(x, xL, xU) = -xabsx_deriv(x) + xL*xabsx_deriv2(x)
+	inflection, flag = newton(xU, 0.0, xU, root_f, root_df, xL, xU)
+	flag && (inflection = golden_section(xL, xU, root_f, xL, xU))
+	if x <= inflection
+		return dline_seg(xabsx, xabsx_deriv, x, xL, inflection)
+	else
+		return xabsx(x), xabsx_deriv(x)
+	end
+end
+function xabsx(x::Interval{Float64})
+	xabsx_xL = Interval(x.lo)*abs(Interval(x.lo))
+	xabsx_xU = Interval(x.hi)*abs(Interval(x.hi))
+	Interval{Float64}(xabsx_xL.lo, xabsx_xU.hi)
+end
+@inline function xabsx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}}
+	xL = x.Intv.lo
+	xU = x.Intv.hi
+	eps_max = xU
+	eps_min = xL
+	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
+	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
+	cc, dcc = cc_xabsx(midcc, xL, xU)
+	cv, dcv = cv_xabsx(midcv, xL, xU)
+	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
+	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
+	cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
+	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
+end
+@inline function xabsx_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Diff}
+	xL = x.Intv.lo
+	xU = x.Intv.hi
+	eps_min = xL
+	eps_max = xU
+	midcc, cc_id = mid3(x.cc, x.cv, eps_max)
+	midcv, cv_id = mid3(x.cc, x.cv, eps_min)
+	cc, dcc = cc_xabsx(midcc, xL, xU)
+	cv, dcv = cv_xabsx(midcv, xL, xU)
+	cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
+	cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
+	return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
+end
+function xabsx(x::MC{N,T}) where {N, T <: RelaxTag}
+	xabsx_kernel(x, xabsx(x.Intv))
+end
+=#
\ No newline at end of file
diff --git a/src/forward_operators/other.jl b/src/forward_operators/other.jl
index 5544da1..260a931 100644
--- a/src/forward_operators/other.jl
+++ b/src/forward_operators/other.jl
@@ -17,6 +17,7 @@ end
 interval_MC(x::MC{S,T}) where {S, T<:RelaxTag} = MC{S,T}(x.Intv)
 
 @inline isnan(x::MC) = isnan(x.cc) || isnan(x.cv)
+@inline isfinite(x::MC) = isfinite(x.cc) && isfinite(x.cv) && isfinite(x.Intv)
 
 ########### Defines differentiable step relaxations
 @inline function cv_step(x::Float64, xL::Float64, xU::Float64)
diff --git a/src/forward_operators/power.jl b/src/forward_operators/power.jl
index 9f210ba..e4993e1 100644
--- a/src/forward_operators/power.jl
+++ b/src/forward_operators/power.jl
@@ -12,10 +12,17 @@
 # defines square operator
 @inline sqr(x::Float64) = x*x
 @inline cv_sqr_NS(x::Float64, xL::Float64, xU::Float64) = x*x
+
 @inline dcv_sqr_NS(x::Float64, xL::Float64, xU::Float64) = 2.0*x
 @inline function cc_sqr(x::Float64, xL::Float64, xU::Float64)
 	if (xU > xL)
-		cc = xL^2 + (xL + xU)*(x - xL)
+		if x == xL
+			cc = xL^2
+		elseif x == xU
+			cc = xU^2
+		else
+			cc = (xL + xU)*x - xU*xL
+		end
 	else
 		cc = xU^2
 	end
diff --git a/src/forward_operators/safe_mode.jl b/src/forward_operators/safe_mode.jl
deleted file mode 100644
index b8c6d27..0000000
--- a/src/forward_operators/safe_mode.jl
+++ /dev/null
@@ -1,78 +0,0 @@
-# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
-# This code is licensed under MIT license (see LICENSE.md for full details)
-#############################################################################
-# McCormick.jl
-# A McCormick operator library in Julia
-# See https://github.com/PSORLab/McCormick.jl
-#############################################################################
-# src/forward_operators/safe_mode.jl
-# UNDER DEVELOPMENT. CURRENTLY NOT USED. Implements a correctly-rounded
-# version of standard McCormick arithmetic.
-#############################################################################
-
-#=
-# Need to add this to Project.toml... requires Julia 1.3+
-using CRLibm
-using RoundingEmulator
-
-struct NSSafe <: RelaxTag end
-
-# Safe addition
-@inline function plus_kernel(x::MC{N,T}, y::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
-	cv = add_down(x.cv, y.cv)
-	cc = add_up(x.cc, y.cc)
-	MC{N,T}(cv, cc, z, x.cv_grad + y.cv_grad, x.cc_grad + y.cc_grad, (x.cnst && y.cnst))
-end
-
-# Safe scalar addition
-@inline function plus_kernel(y::Float64, x::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
-	cv = add_down(x.cv, y)
-	cc = add_up(x.cc, y)
-	MC{N,T}(cv, cc, z, x.cv_grad, x.cc_grad, x.cnst)
-end
-
-# Safe subtraction
-@inline function minus_kernel(x::MC{N,T}, y::MC{N,T}, z::Interval{Float64}) where {N, T <:  Union{NSSafe}}
-	cv = sub_down(x.cv, y.cc)
-	cc = sub_up(x.cc, y.cv)
-	MC{N,T}(cv, cc, z, x.cv_grad - y.cc_grad, x.cc_grad - y.cv_grad, (x.cnst && y.cnst))
-end
-
-# Safe scalar subtraction
-@inline function minus_kernel(x::MC{N,T}, c::Float64, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
-	cv = sub_down(x.cv, c)
-	cc = sub_up(x.cc, c)
-	MC{N,T}(cv, cc, z, x.cv_grad, x.cc_grad, x.cnst)
-end
-@inline function minus_kernel(c::Float64, x::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
-	cv = sub_down(c, x.cc)
-	cc = sub_up(c, x.cv)
-	MC{N,T}(cv, cc, z, -x.cc_grad, -x.cv_grad, x.cnst)
-end
-
-# Safe Scalar Multiplication
-@inline function mult_kernel(x::MC{N,T}, c::Float64, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
-	if c >= 0.0
-		cv = mul_down(c, x.cv)
-		cc = mul_up(c, x.cc)
-		zMC = MC{N,T}(cv, cc, z, c*x.cv_grad, c*x.cc_grad, x.cnst)
-	else
-		cv = mul_down(c, x.cc)
-		cc = mul_up(c, x.cv)
-		zMC = MC{N,T}(cc, cv, z, c*x.cc_grad, c*x.cv_grad, x.cnst)
-	end
-	return zMC
-end
-
-# Scalar Safe Kernel
-@inline function div_kernel(x::MC{N,T}, c::Float64, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
-	if c >= 0.0
-		cv = div_down(x.cv, c)
-		cc = div_up(x.cc, c)
-	else
-		cc = div_down(x.cv, c)
-		cv = div_up(x.cc, c)
-	end
-	return MC{N,T}(cc, cv, z, c*x.cc_grad, c*x.cv_grad, x.cnst)
-end
-=#
diff --git a/src/forward_operators/trilinear.jl b/src/forward_operators/trilinear.jl
index 4dae025..c4e78bc 100644
--- a/src/forward_operators/trilinear.jl
+++ b/src/forward_operators/trilinear.jl
@@ -1437,6 +1437,7 @@ end
     z = x1.Intv*x2.Intv*x3.Intv
 	return trilinear_kernel(x1, x2, x3, z)
 end
+trilinear(x, y, z) = x*y*z
 
 #=
 *(a::MC{N,T}, b::MC{N,T}, c::MC{N,T}, d::MC{N,T}) where {N, T<:Union{NS,MV}} = *(*(a,b,c),d)
diff --git a/test/forward_mccormick.jl b/test/forward_mccormick.jl
index 25499a9..16429c8 100644
--- a/test/forward_mccormick.jl
+++ b/test/forward_mccormick.jl
@@ -208,35 +208,35 @@ if ~(VERSION < v"1.1-")
       y12 = gelu(x12)
 
       @test isapprox(gelu(7), 6.999999999991041, atol=1E-8)
-      @test isapprox(gelu(-0.75), -0.16997051428265114, atol=1E-8)
+      @test isapprox(gelu(-0.75), -0.17003944483437966, atol=1E-8)
 
-      @test isapprox(y1.cv, 0.053982783727702904, atol=1E-8)
-      @test isapprox(y1.cc, 0.6997891980967129, atol=1E-8)
-      @test isapprox(y2.cv, -0.15426876936299344, atol=1E-8)
-      @test isapprox(y2.cc, 0.3997891980967129, atol=1E-8)
-      @test isapprox(y3.cv, 0.5306254434438489, atol=1E-8)
-      @test isapprox(y3.cc, 0.5565428241289478, atol=1E-8)
+      @test isapprox(y1.cv, 0.053982751045435165, atol=1E-8)
+      @test isapprox(y1.cc, 0.6995715769802328, atol=1E-8)
+      @test isapprox(y2.cv, -0.15428599017485609, atol=1E-8)
+      @test isapprox(y2.cc, 0.3995715769802329, atol=1E-8)
+      @test isapprox(y3.cv, 0.5305701347051167, atol=1E-8)
+      @test isapprox(y3.cc, 0.5564855232561616, atol=1E-8)
 
-      @test isapprox(y4.cv, -5.222141651051171e-8, atol=1E-8)
+      @test isapprox(y4.cv, -2.9638201337611747e-9, atol=1E-8)
       @test isapprox(y4.cc, -2.610399119085116e-10, atol=1E-8)
-      @test isapprox(y5.cv, 6.4999999477785835, atol=1E-8)
+      @test isapprox(y5.cv, 6.49999999703618, atol=1E-8)
       @test isapprox(y5.cc, 6.49999999973896, atol=1E-8)
-      @test isapprox(y6.cv, 0.5306254434438489, atol=1E-8)
-      @test isapprox(y6.cc, 0.5565428241289478, atol=1E-8)
-
-      @test isapprox(y7.cv, 3.073809840032428, atol=1E-8)
-      @test isapprox(y7.cc, 3.299185798223376, atol=1E-8)
-      @test isapprox(y8.cv, 0.9505770065222403, atol=1E-8)
-      @test isapprox(y8.cc, 1.697372160948963, atol=1E-8)
-      @test isapprox(y9.cv, 0.18537342665668577, atol=1E-8)
-      @test isapprox(y9.cc, 1.4368229984446212, atol=1E-8)
-
-      @test isapprox(y10.cv, -0.14923266704102095, atol=1E-8)
-      @test isapprox(y10.cc, -0.02584062996129348, atol=1E-8)
-      @test isapprox(y11.cv, 0.7343458871879165, atol=1E-8)
-      @test isapprox(y11.cc, 0.8328461064510111, atol=1E-8)
-      @test isapprox(y12.cv, 0.18537342665668577, atol=1E-8)
-      @test isapprox(y12.cc, 0.5343290419885047, atol=1E-8)
+      @test isapprox(y6.cv, 0.5305701347051167, atol=1E-8)
+      @test isapprox(y6.cc, 0.5564855232561616, atol=1E-8)
+
+      @test isapprox(y7.cv, 3.0739573701225944, atol=1E-8)
+      @test isapprox(y7.cc, 3.299383802344626, atol=1E-8)
+      @test isapprox(y8.cv, 0.9504228742199268, atol=1E-8)
+      @test isapprox(y8.cc, 1.6977056770764198, atol=1E-8)
+      @test isapprox(y9.cv, 0.18537092354275922, atol=1E-8)
+      @test isapprox(y9.cc, 1.4369627939042764, atol=1E-8)
+
+      @test isapprox(y10.cv, -0.14941895382548911, atol=1E-8)
+      @test isapprox(y10.cc, -0.026070987682690627, atol=1E-8)
+      @test isapprox(y11.cv, 0.7342284887097207, atol=1E-8)
+      @test isapprox(y11.cc, 0.8326705877810306, atol=1E-8)
+      @test isapprox(y12.cv, 0.18537092354275922, atol=1E-8)
+      @test isapprox(y12.cc, 0.534168383326423, atol=1E-8)
    end
 end
 
@@ -689,15 +689,15 @@ end
    @test isnan(Y1^(-3.0))
    @test isnan(Y2^(-3.0))
 
-   @test isapprox(McCormick.tanh_deriv(3.0, 2.0, 4.0), 0.009866, rtol=1E-3)
+   @test isapprox(McCormick.tanh_deriv(3.0), 0.009866, rtol=1E-3)
    @test isapprox(McCormick.tanh_envd(3.0, 2.0, 4.0), 6.258589, rtol=1E-3)
-   @test isapprox(McCormick.atan_deriv(3.0, 2.0, 4.0), 0.100, rtol=1E-3)
-   @test isapprox(McCormick.asin_deriv(0.5, 0.0, 0.9), 1.1547005, rtol=1E-3)
-   @test isapprox(McCormick.tan_deriv(3.0, 2.0, 4.0), 1.0203195, rtol=1E-3)
+   @test isapprox(McCormick.atan_deriv(3.0), 0.100, rtol=1E-3)
+   @test isapprox(McCormick.asin_deriv(0.5), 1.1547005, rtol=1E-3)
+   @test isapprox(McCormick.tan_deriv(3.0), 1.0203195, rtol=1E-3)
    @test isapprox(McCormick.tan_envd(3.0, 2.0, 4.0), -0.570704, rtol=1E-3)
-   @test isapprox(McCormick.acos_deriv(0.5, 0.0, 0.9), -1.1547005, rtol=1E-3)
+   @test isapprox(McCormick.acos_deriv(0.5), -1.1547005, rtol=1E-3)
    @test isapprox(McCormick.acos_env(0.5, 0.0, 0.9), -0.04655, rtol=1E-3)
-   @test isapprox(McCormick.asinh_deriv(3.0, 2.0, 4.0), 0.31622, rtol=1E-3)
+   @test isapprox(McCormick.asinh_deriv(3.0), 0.31622, rtol=1E-3)
 
    xD = MC{2,Diff}(2.0,2.0,Interval{Float64}(1.0,4.0), seed_gradient(2,Val(2)), seed_gradient(2,Val(2)),false)
    xNS = MC{2,NS}(2.0,2.0,Interval{Float64}(1.0,4.0), seed_gradient(2,Val(2)), seed_gradient(2,Val(2)),false)
@@ -1807,7 +1807,7 @@ end
 
    @test isapprox(sigmoid(7), 0.9990889488055994, atol=1E-8)
    @test isapprox(sigmoid(0.75), 0.679178699175393, atol=1E-8)
-   @test isapprox(McCormick.sigmoid_deriv2(0.75), 0.4003874719875089, atol=1E-8)
+   @test isapprox(McCormick.sigmoid_deriv2(0.75), -0.07808428307814444, atol=1E-8)
 
    @test isapprox(y1.cv, 0.5084495314200309, atol=1E-8)
    @test isapprox(y1.cc, 0.5357111749761996, atol=1E-8)
@@ -1859,24 +1859,24 @@ end
    x11 = MC{1,NS}(0.9, Interval(-0.8, 1.1), 1)
    x12 =  MC{1,NS}(0.3, Interval(-0.8, 1.3), 1)
 
-   y1 = swish1(x1)
-   y2 = swish1(x2)
-   y3 = swish1(x3)
+   y1 = swish(x1)
+   y2 = swish(x2)
+   y3 = swish(x3)
 
-   y4 = swish1(x4)
-   y5 = swish1(x5)
-   y6 = swish1(x6)
+   y4 = swish(x4)
+   y5 = swish(x5)
+   y6 = swish(x6)
 
-   y7 = swish1(x7)
-   y8 = swish1(x8)
-   y9 = swish1(x9)
+   y7 = swish(x7)
+   y8 = swish(x8)
+   y9 = swish(x9)
 
-   y10 = swish1(x10)
-   y11 = swish1(x11)
-   y12 = swish1(x12)
+   y10 = swish(x10)
+   y11 = swish(x11)
+   y12 = swish(x12)
 
-   @test isapprox(swish1(7), 6.993622641639196, atol=1E-8)
-   @test isapprox(swish1(-0.75), -0.24061597561845527, atol=1E-8)
+   @test isapprox(swish(7), 6.993622641639196, atol=1E-8)
+   @test isapprox(swish(-0.75), -0.24061597561845527, atol=1E-8)
 
    @test isapprox(y1.cv, 0.052497918747894, atol=1E-8)
    @test isapprox(y1.cc, 0.5263617142904655, atol=1E-8)
@@ -1906,7 +1906,7 @@ end
    @test isapprox(y12.cv, 0.1723327550434977, atol=1E-8)
    @test isapprox(y12.cc, 0.4170112431680709, atol=1E-8)
 
-   @test isapprox(McCormick.swish1_env(1.0, 2.0, 3.0), -0.1028650654542731, atol=1E-8)
+   @test isapprox(McCormick.swish_env(1.0, 2.0, 3.0), -0.1028650654542731, atol=1E-8)
 
    x4 = MC{1,NS}(-6.5, Interval(-7.5, -5.5), 1)
    x5 = MC{1,NS}(6.5, Interval(5.5, 7.5), 1)