Refactor, remove artificial restrictions

src/mittag_leffler.jl

@@ -4,7 +4,7 @@
 
 # When support for v0.6 is dropped, use fully qualified name
 if VERSION >= v"0.7-"
-    import SpecialFunctions: gamma
+    import SpecialFunctions: gamma, erfc
 end
 
 macro br(n)
@@ -91,7 +91,7 @@ function mittleffsum2(α,β,z,ρ)
     @br 2
     k0 = max(ceil(Int,(1-β)/α), ceil(Int, log(ρ*(1-abs(z)))/log(abs(z))))
     s = zero(z)
-    for k=0:k0
+    for k=0:(k0)
         s += z^k/gamma(β+α*k)
     end
     return s
@@ -144,63 +144,91 @@ Compute the Mittag-Leffler function at `z` for parameters `α,β` with
 accuracy `ρ`.
 """
 function mittlefferr(α,β,z,ρ)
-    ρ > 0 || throw(DomainError())
+    ρ > 0 || throw(DomainError(ρ))
     _mittlefferr(α,β,z,ρ)
 end
 
-_mittlefferr(α::Real,β::Real,z::Real,ρ::Real) = real(_mittleff(α,β,z,ρ))
-_mittlefferr(α::Real,β::Real,z::Complex,ρ::Real) = _mittleff(α,β,z,ρ)
+_mittlefferr(α::Real, β::Real, z::Real, ρ::Real) = real(_mittleff(α, β, z, ρ))
+_mittlefferr(α::Real, β::Real, z::Complex, ρ::Real) = _mittleff(α, β, z, ρ)
 
 # The second definition would work for both complex and real
 myeps(x) = x |> one |> float |> eps
 myeps(x::Complex) =  x |> real |> myeps
 
-
 """
     mittleff(α,β,z)
 
 Compute the Mittag-Leffler function at `z` for parameters `α,β`.
 """
-mittleff(α, β, z) = _mittlefferr(α,β,z,myeps(z))
+mittleff(α, β, z) = _mittleff(α,β,z)
+#mittleff(α, β, z) = _mittlefferr(α,β,z,myeps(z))
 mittleff(α, β, z::Union{Integer,Complex{T}}) where {T<:Integer} = mittleff(α, β, float(z))
 
 """
     mittleff(α,z)
 
 Compute `mittleff(α,1,z)`.
 """
-mittleff(α,z) = _mittlefferr(α,1,z,myeps(z))
-mittleff(α,z::Union{Integer,Complex{T}}) where {T<:Integer} = mittleff(α, float(z))
+mittleff(α, z) = _mittleff(α,1,z)
+mittleff(α, z::Union{Integer,Complex{T}}) where {T<:Integer} = mittleff(α, float(z))
+
+function _mittleff_special_beta_one(α,z)
+    z == 0 && return one(z)
+    α == 1/2 && return exp(z^2)*erfc(-z)
+    α == 0 && return 1/(1-z)  # FIXME needs domain check
+    α == 1 && return exp(z)
+    zc = isreal(z) && z < 0 ? complex(z) : z # Julia sometimes requires explicit complex numbers for efficiency
+    α == 2 && return cosh(sqrt(zc))
+    α == 3 && return (1//3)*(exp(zc^(1//3)) + 2*exp(-zc^(1//3)/2) * cos(sqrt(convert(typeof(zc),3))/2 * zc^(1//3)))
+    α == 4 && return (1//2)*(cosh(zc^(1//4)) + cos(zc^(1//4)))
+    return nothing
+end
 
-function _mittleff(α,β,z,ρ)
-#    if β == 1
-        # if α == 1/2   # disable this. There is never an error here. But, this triggers a mysterious, untraceable bug in quadgk
-        #     res = try
-        #         exp(z^2)*erfc(-z)
-        #     catch
-        #         error("Failed in exp. erfc")
-        #     end
-        #     return res
-        # end
-        # Disable these, because we would have to take care of domain errors, etc.
-        # α == 0 && return 1/(1-z)
-        # α == 1 && return exp(z)
-        # α == 2 && return cosh(sqrt(z))
-        # α == 3 && return (1//3)*(exp(z^(1//3)) + 2*exp(-z^(1//3)/2) * cos(sqrt(convert(typeof(z),3))/2 * z^(1//3)))
-        # α == 4 && return (1//2)*(cosh(z^(1//4)) + cos(z^(1//4)))
-#    end
+# This tries to use small unions in an efficient way, but does not.
+# It adds at least 10 ns
+# So, we do not use it at the moment.
+function _mittleff_no_eps(α, β, z)
     z == 0 && return 1/gamma(β)
-    α == 1 && β == 1 && return(exp(z))
-    α < 0  && throw(DomainError())
-    az = abs(z)
     1 < α && return mittleffsum(α,β,z)
+    return nothing
+end
+
+function _mittleff_slow_with_eps(α, β, z, ρ)
+    az = abs(z)
     az < 1 && return mittleffsum2(α,β,z,ρ)
     az > floor(10+5*α) && return choosesum(α,β,z,ρ)
     return mittleffints(α,β,z,ρ)
 end
 
-_mittleff(α,β,z) = mittleff(α,β,z,myeps(z))
+_mittleff(α::Real, β::Real, z::Real) = real(_mittleff0(α, β, z))
+_mittleff(α, β, z) = _mittleff0(α, β, z)
+
+function _mittleff0(α, β, z)
+    if β == 1
+        res = _mittleff_special_beta_one(α,z)
+        res != nothing && return res
+    end
+    z == 0 && return 1/gamma(β)
+    if β == 2
+        α == 1 && return (exp(z) - 1)/z
+        zc = z < 0 ? complex(z) : z
+        α == 2 && return sinh(sqrt(zc))/sqrt(zc)
+    end
+    1 < α && return mittleffsum(α,β,z)
+    ρ = myeps(z)
+    return _mittleff_slow_with_eps(α, β, z, ρ)
+end
 
+function _mittleff(α,β,z,ρ)
+    if β == 1
+        res = _mittleff_special_beta_one(α,z)
+        res != nothing && return res
+    end
+    z == 0 && return 1/gamma(β)
+    α == 1 && β == 1 && return(exp(z))
+    1 < α && return mittleffsum(α,β,z)
+    return _mittleff_slow_with_eps(α, β, z, ρ)
+end
 
 """
     mittleffderiv(α,β,z)