diff --git a/Project.toml b/Project.toml
index f8a6b7f..1bd08a4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,4 +23,4 @@ ForwardDiff = "0.10"
 Setfield = "1.1"
 StaticArrays = "1.6"
 julia = "1.8"
-Enzyme = "0.11.20"
+Enzyme = "0.11"
diff --git a/README.md b/README.md
index 0166dcf..eeda8c8 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,5 @@
 # JacobiElliptic
-JacobiElliptic is an implementation of Toshio Fukushima's algorithms for calculating [Elliptic Integrals and Jacobi Elliptic Functions](https://ieeexplore.ieee.org/document/7203795). 
+JacobiElliptic is an implementation of Toshio Fukushima's & Billie C. Carlson's for calculating [Elliptic Integrals and Jacobi Elliptic Functions](https://ieeexplore.ieee.org/document/7203795). 
 
 ## Features
   - Type stable and preserving
diff --git a/ext/JacobiEllipticEnzymeExt.jl b/ext/JacobiEllipticEnzymeExt.jl
index 892edd9..405ebc9 100644
--- a/ext/JacobiEllipticEnzymeExt.jl
+++ b/ext/JacobiEllipticEnzymeExt.jl
@@ -18,72 +18,61 @@ function ∂F_∂ϕ(ϕ, m)
 end
 
 function forward(
-    ::Const{typeof(JacobiElliptic.F)},
-    ::Type,
+    func::Const{typeof(JacobiElliptic.F)},
+    ::Type{<:Duplicated},
     ϕ::Const,
     m::Duplicated
 ) 
-    println("In custom forward rule.")
-
-    return ∂F_∂m(ϕ.val, m.val)*m.dval
+    return Duplicated(func.val(ϕ.val, m.val), ∂F_∂m(ϕ.val, m.val)*m.dval)
 end
 
 function forward(
-    ::Const{typeof(JacobiElliptic.F)},
-    ::Type,
+    func::Const{typeof(JacobiElliptic.F)},
+    ::Type{<:Duplicated},
     ϕ::Duplicated,
     m::Const
 ) 
-    println("In custom forward rule.")
-
-    return ∂F_∂ϕ(ϕ.val, m.val)*ϕ.dval
+    return Duplicated(func.val(ϕ.val, m.val), ∂F_∂ϕ(ϕ.val, m.val)*ϕ.dval)
 end
 
 function forward(
-    ::Const{typeof(JacobiElliptic.F)},
-    ::Type,
+    func::Const{typeof(JacobiElliptic.F)},
+    ::Type{<:Duplicated},
     ϕ::Duplicated,
     m::Duplicated
 ) 
-    println("In custom forward rule.")
-
-    return ∂F_∂m(ϕ.val, m.val)*m.dval + ∂F_∂ϕ(ϕ.val, m.val)*ϕ.dval
+    return Duplicated(func.val(ϕ.val, m.val), ∂F_∂m(ϕ.val, m.val)*m.dval + ∂F_∂ϕ(ϕ.val, m.val)*ϕ.dval)
 end
 
 function augmented_primal(
-    config,
+    config::ConfigWidth{N},
     func::Const{typeof(JacobiElliptic.F)},
     ::Type{<:Active},
-    tape,
-    ϕ,
-    m
- ) 
-    println("In custom augmented primal rule.")
+    ϕ::Union{Const,Active},
+    m::Union{Const,Active}
+ ) where {N}
+    #println("In custom augmented primal rule.")
     # Save x in tape if x will be overwritten
-    primal = EnzymeRules.needs_primal(config) ? func(ϕ.val, m.val) : nothing
+    primal = EnzymeRules.needs_primal(config) ? func.val(ϕ.val, m.val) : nothing
 
-    return EnzymeRules.AugmentedReturn(primal, nothing, tape)
+    return EnzymeRules.AugmentedReturn(primal, nothing, nothing)
 end
 
-
-function reverse(config::ConfigWidth{N}, ::Const{JacobiElliptic.F}, dret::Active, tape, ϕ::Const, m::Active) where {N}  
-    println("In custom reverse rule.")
+function reverse(config::ConfigWidth{1}, ::Const{typeof(JacobiElliptic.F)}, dret::Active, tape, ϕ::Const, m::Active) 
     # retrieve x value, either from original x or from tape if x may have been overwritten.
     mval = m.val
     dm = ∂F_∂m(ϕ.val, mval) * dret.val
-    return (dm, )
+    return (nothing, dm)
 end
 
-function reverse(config::ConfigWidth{N}, ::Const{JacobiElliptic.F}, dret::Active, tape, ϕ::Active, m::Const) where {N}  
-    println("In custom reverse rule.")
+function reverse(config::ConfigWidth{1}, ::Const{typeof(JacobiElliptic.F)}, dret::Active, tape, ϕ::Active, m::Const) 
     # retrieve x value, either from original x or from tape if x may have been overwritten.
-    ϕval = EnzymeRules.overwritten(config)[2] ? tape : ϕ.val
+    ϕval = ϕ.val
     dϕ = ∂F_∂ϕ(ϕval, m.val) * dret.val
-    return (dϕ, )
+    return (dϕ, nothing)
 end
 
-function reverse(config::ConfigWidth{N}, ::Const{JacobiElliptic.F}, dret::Active, tape, ϕ::Active, m::Active) where {N}  
-    println("In custom reverse rule.")
+function reverse(config, ::Const{typeof(JacobiElliptic.F)}, dret::Active, tape, ϕ::Union{Active, Duplicated}, m::Union{Active, Duplicated})
     # retrieve x value, either from original x or from tape if x may have been overwritten.
     ϕval = ϕ.val
     mval = m.val
@@ -92,5 +81,4 @@ function reverse(config::ConfigWidth{N}, ::Const{JacobiElliptic.F}, dret::Active
     return (dϕ, dm)
 end
 
-
 end # module
\ No newline at end of file
diff --git a/src/Fukushima.jl b/src/Fukushima.jl
index 22915e4..7a5b6c1 100644
--- a/src/Fukushima.jl
+++ b/src/Fukushima.jl
@@ -19,7 +19,6 @@ Returns the complete elliptic integral of the first kind.
 - `m` : Elliptic modulus
 """
 function K(m::T) where T
-    !(-1 ≤ m ≤ 1) && throw(DomainError("argument m not in [0,1]"))
     m == one(T) && return T(Inf)
 
     # I didn't really see any speedup from evalpoly, so I left the evaluation in this form
@@ -50,7 +49,7 @@ function K(m::T) where T
     x == zero(T) && return T(π/2)
     x == one(T) && return T(Inf)
     x > one(T) && return T(NaN)
-    
+        
     if x < T(0.1)
         t = poly1( x - T(0.05));
     elseif ( x < T(0.2)) 
@@ -79,7 +78,7 @@ function K(m::T) where T
     end
     # Complete the transformation mentioned above for m < 0:
     flag && return t / sqrt( one(T) - m );
-    
+        
     return t
 end
 
@@ -93,7 +92,6 @@ Returns the complete elliptic integral of the second kind.
 - `m` : Elliptic modulus
 """
 function E(m::T) where T
-    !(-1 ≤ m ≤ 1) && throw(DomainError("argument m not in [0,1]"))
     m == zero(T) && return T(π/2)
     m == one(T) && return one(T)
 
@@ -224,13 +222,14 @@ function asn(s::A, m::B) where {A,B}
     yA = T(0.04094) - T(0.00652)*m
     y = s * s
     y < yA && return s*serf(y, m)
-
+    
     p = one(T)
     for _ in 1:10
         y = y * inv((1+√(1-y))*(1+√(1-m*y)))
         p += p
         y < yA && return p*√y*serf(y, m)
     end
+    return T(NaN)
 end
 
 """
@@ -257,6 +256,7 @@ function acn(c::A, m::B) where {A,B}
         x = (√x + d)/(1+d)
         p += p
     end
+    return T(NaN)
 end
 
 @inline function rawF(sinφ::A, m::B) where {A,B}
@@ -552,7 +552,7 @@ function Pi(n::A, φ::B, m::C) where {A,B,C}
         return (FukushimaT(t1, h1) - n1*J(n1, φ, m))
     end
     return n*J(n, φ, m) + F(φ, m)
-end
+    end
 
 """
 ``J (n;\\varphi \\,|\\,m)=\\frac{\\Pi(n;\\varphi|\\, m) - F(\\varphi|\\,m)}{n}.``
@@ -566,7 +566,7 @@ Returns the associate incomplete elliptic integral of the third kind.
 - `m` : Elliptic modulus
 """
 function J(n::A, φ::B, m::C) where {A, B, C} #Appendix A
-    T = promote_type(A,B,C)
+T = promote_type(A,B,C)
     # Reduction of Amplitude
     φ == zero(T) && return zero(T)
     φ == T(π/2) && return J(n,m)
@@ -818,16 +818,16 @@ function FukushimaT(t::A, h::B) where {A,B}
 		return t
 	else
         arg = t * √(-h)
-        ans = abs(arg) < one(T) ? atanh(arg) : custom_atanh(arg)
-        return ans / √(-h)
-	end
+                    ans = abs(arg) < one(T) ? atanh(arg) : custom_atanh(arg)
+            return ans / √(-h)
+        	end
 end
 
 #https://link-springer-com.ezp-prod1.hul.harvard.edu/article/T(10).1007/BF02165405
 function cel(kc::A, p::B, a::C, b::D) where {A,B,C,D}
     T = promote_type(A,B,C,D)
     #ca = T(1e-6)
-    ca = eps(T)
+ca = eps(T)
     kc = abs(kc)
     e = kc
     m = one(T)
@@ -932,7 +932,7 @@ function fold_1_00(u1::T, m::T, Kscreen::T, Kactual::T, kp::T) where T
 
     if u1 > Kscreen/2
         sn, cn, dn = fold_0_50(Kactual - u1, m, Kscreen, Kactual, kp)
-    elseif u1 > Kscreen/4
+    elseif u1 > Kscreen/4 
         sn, cn, dn = fold_0_25(Kactual/2 - u1, m, kp)
     else
         sn, cn, dn = _ΔXNloop(u1, m, u1 > zero(T) ?  max(6+(floor(log2(u1))), one(T)) : zero(T))
diff --git a/src/JacobiElliptic.jl b/src/JacobiElliptic.jl
index 26ea2f9..8435c11 100644
--- a/src/JacobiElliptic.jl
+++ b/src/JacobiElliptic.jl
@@ -22,7 +22,7 @@ struct Fukushima <: AbstractAlgorithm end
 struct Carlson <: AbstractAlgorithm end
 
 
-func_syms = [:E, :F, :K, :Pi, :J, :sn, :cn, :dn, :nn, :sd, :dd, :nd, :sc, :cc, :dc, :nc, :ss, :cs, :ds, :ns, :am, :cd]
+func_syms = [:E, :F, :K, :Pi, :J, :sn, :cn, :dn, :nn, :sd, :dd, :nd, :sc, :cc, :dc, :nc, :ss, :cs, :ds, :ns, :cd]
 sym_list = []
 
 
@@ -46,5 +46,6 @@ end
 
 asn = FukushimaAlg.asn
 acn = FukushimaAlg.asn
+am = CarlsonAlg.am
 
 end
\ No newline at end of file