improve support for custom Number types (#50)

Author: aplavin

Project.toml

@@ -1,6 +1,6 @@
 name = "InverseFunctions"
 uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.14"
+version = "0.1.15"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -18,6 +18,7 @@ julia = "1"
 [extras]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["Dates", "Documenter"]
+test = ["Dates", "Documenter", "Unitful"]

src/functions.jl

@@ -5,57 +5,46 @@
 
 Inverse of `sqrt(x)` for non-negative `x`.
 """
-square(x) = x^2
-function square(x::Real)
-    x < zero(x) && throw(DomainError(x, "`square` is defined as the inverse of `sqrt` and can only be evaluated for non-negative values"))
+function square(x)
+    if is_real_type(typeof(x)) && x < zero(x)
+        throw(DomainError(x, "`square` is defined as the inverse of `sqrt` and can only be evaluated for non-negative values"))
+    end
     return x^2
 end
 
 
-function invpow2(x::Real, p::Integer)
-    if x ≥ zero(x) || isodd(p)
-        copysign(abs(x)^inv(p), x)
+function invpow_arg2(x::Number, p::Real)
+    if is_real_type(typeof(x))
+        x ≥ zero(x) ? x^inv(p) :  # x > 0 - trivially invertible
+            isinteger(p) && isodd(Integer(p)) ? copysign(abs(x)^inv(p), x) :  # p odd - invertible even for x < 0
+            throw(DomainError(x, "inverse for x^$p is not defined at $x"))
     else
-        throw(DomainError(x, "inverse for x^$p is not defined at $x"))
-    end
-end
-function invpow2(x::Real, p::Real)
-    if x ≥ zero(x)
-        x^inv(p)
-    else
-        throw(DomainError(x, "inverse for x^$p is not defined at $x"))
-    end
-end
-function invpow2(x, p::Real)
-    # complex x^p is only invertible for p = 1/n
-    if isinteger(inv(p))
-        x^inv(p)
-    else
-        throw(DomainError(x, "inverse for x^$p is not defined at $x"))
+        # complex x^p is invertible only for p = 1/n
+        isinteger(inv(p)) ? x^inv(p) : throw(DomainError(x, "inverse for x^$p is not defined at $x"))
     end
 end
 
-function invpow1(b::Real, x::Real)
+function invpow_arg1(b::Real, x::Real)
     if b ≥ zero(b) && x ≥ zero(x)
         log(b, x)
     else
         throw(DomainError(x, "inverse for $b^x is not defined at $x"))
     end
 end
 
-function invlog1(b::Real, x::Real)
+function invlog_arg1(b::Real, x::Real)
     if b ≥ zero(b)
         b^x
     else
         throw(DomainError(x, "inverse for log($b, x) is not defined at $x"))
     end
 end
-invlog1(b, x) = b^x
+invlog_arg1(b::Number, x::Number) = b^x
 
-invlog2(b, x) = x^inv(b)
+invlog_arg2(b::Number, x::Number) = x^inv(b)
 
 
-function invdivrem((q, r), divisor)
+function invdivrem((q, r)::NTuple{2,Number}, divisor::Number)
     res = muladd(q, divisor, r)
     if abs(r) ≤ abs(divisor) && (iszero(r) || sign(r) == sign(res))
         res
@@ -64,10 +53,18 @@ function invdivrem((q, r), divisor)
     end
 end
 
-function invfldmod((q, r), divisor)
+function invfldmod((q, r)::NTuple{2,Number}, divisor::Number)
     if abs(r) ≤ abs(divisor) && (iszero(r) || sign(r) == sign(divisor))
         muladd(q, divisor, r)
     else
         throw(DomainError((q, r), "inverse for fldmod(x) is not defined at this point"))
     end
 end
+
+
+# check if T is a real-Number type
+# this is not the same as T <: Real which immediately excludes custom Number subtypes such as unitful numbers
+# also, isreal(x) != is_real_type(typeof(x)): the former is true for complex numbers with zero imaginary part
+@inline is_real_type(@nospecialize _::Type{<:Real}) = true
+@inline is_real_type(::Type{T}) where {T<:Number} = real(T) == T
+@inline is_real_type(@nospecialize _::Type) = false

src/inverse.jl

@@ -159,14 +159,14 @@ inverse(::typeof(sqrt)) = square
 inverse(::typeof(square)) = sqrt
 
 inverse(::typeof(cbrt)) = Base.Fix2(^, 3)
-inverse(f::Base.Fix2{typeof(^)}) = iszero(f.x) ? throw(DomainError(f.x, "Cannot invert x^$(f.x)")) : Base.Fix2(invpow2, f.x)
-inverse(f::Base.Fix2{typeof(invpow2)}) = Base.Fix2(^, f.x)
-inverse(f::Base.Fix1{typeof(^)}) = Base.Fix1(invpow1, f.x)
-inverse(f::Base.Fix1{typeof(invpow1)}) = Base.Fix1(^, f.x)
-inverse(f::Base.Fix1{typeof(log)}) = Base.Fix1(invlog1, f.x)
-inverse(f::Base.Fix1{typeof(invlog1)}) = Base.Fix1(log, f.x)
-inverse(f::Base.Fix2{typeof(log)}) = Base.Fix2(invlog2, f.x)
-inverse(f::Base.Fix2{typeof(invlog2)}) = Base.Fix2(log, f.x)
+inverse(f::Base.Fix2{typeof(^)}) = iszero(f.x) ? throw(DomainError(f.x, "Cannot invert x^$(f.x)")) : Base.Fix2(invpow_arg2, f.x)
+inverse(f::Base.Fix2{typeof(invpow_arg2)}) = Base.Fix2(^, f.x)
+inverse(f::Base.Fix1{typeof(^)}) = Base.Fix1(invpow_arg1, f.x)
+inverse(f::Base.Fix1{typeof(invpow_arg1)}) = Base.Fix1(^, f.x)
+inverse(f::Base.Fix1{typeof(log)}) = Base.Fix1(invlog_arg1, f.x)
+inverse(f::Base.Fix1{typeof(invlog_arg1)}) = Base.Fix1(log, f.x)
+inverse(f::Base.Fix2{typeof(log)}) = Base.Fix2(invlog_arg2, f.x)
+inverse(f::Base.Fix2{typeof(invlog_arg2)}) = Base.Fix2(log, f.x)
 
 inverse(f::Base.Fix2{typeof(divrem)}) = Base.Fix2(invdivrem, f.x)
 inverse(f::Base.Fix2{typeof(invdivrem)}) = Base.Fix2(divrem, f.x)

test/test_inverse.jl

@@ -2,6 +2,7 @@
 
 using Test
 using InverseFunctions
+using Unitful
 using Dates
 
 
@@ -79,14 +80,19 @@ end
     end
 
     # ensure that inverses have domains compatible with original functions
+    @test_throws DomainError inverse(sqrt)(-1.0)
+    InverseFunctions.test_inverse(sqrt, complex(-1.0))
+    InverseFunctions.test_inverse(sqrt, complex(1.0))
     @test_throws DomainError inverse(Base.Fix1(*, 0))
     @test_throws DomainError inverse(Base.Fix2(^, 0))
     @test_throws DomainError inverse(Base.Fix1(log, -2))(5)
     @test_throws DomainError inverse(Base.Fix1(log, 2))(-5)
     InverseFunctions.test_inverse(inverse(Base.Fix1(log, 2)), complex(-5))
     @test_throws DomainError inverse(Base.Fix2(^, 0.5))(-5)
     @test_throws DomainError inverse(Base.Fix2(^, 0.51))(complex(-5))
+    @test_throws DomainError inverse(Base.Fix2(^, 2))(complex(-5))
     InverseFunctions.test_inverse(Base.Fix2(^, 0.5), complex(-5))
+    InverseFunctions.test_inverse(Base.Fix2(^, -1), complex(-5.))
     @test_throws DomainError inverse(Base.Fix2(^, 2))(-5)
     @test_throws DomainError inverse(Base.Fix1(^, 2))(-5)
     @test_throws DomainError inverse(Base.Fix1(^, -2))(3)
@@ -130,6 +136,20 @@ end
     end
 end
 
+@testset "unitful" begin
+    # the majority of inverse just propagate to underlying mathematical functions and don't have any issues with unitful numbers
+    # only those that behave treat real numbers differently have to be tested here
+    x = rand()u"m"
+    InverseFunctions.test_inverse(sqrt, x)
+    @test_throws DomainError inverse(sqrt)(-x)
+
+    InverseFunctions.test_inverse(Base.Fix2(^, 2), x)
+    @test_throws DomainError inverse(Base.Fix2(^, 2))(-x)
+    InverseFunctions.test_inverse(Base.Fix2(^, 3), x)
+    InverseFunctions.test_inverse(Base.Fix2(^, 3), -x)
+    InverseFunctions.test_inverse(Base.Fix2(^, -3.5), x)
+end
+
 @testset "dates" begin
     InverseFunctions.test_inverse(Dates.date2epochdays, Date(2020, 1, 2); compare = ===)
     InverseFunctions.test_inverse(Dates.datetime2epochms, DateTime(2020, 1, 2, 12, 34, 56); compare = ===)