fix inverse(^) domain

src/functions.jl

@@ -13,29 +13,22 @@ end
 
 
 function invpow2(x::Real, p::Integer)
-    if x ≥ zero(x) || isodd(p)
-        copysign(abs(x)^inv(p), x)
-    else
-        throw(DomainError(x, "inverse for x^$p is not defined at $x"))
-    end
+    # exception should never happen in actual use: this check is done in inverse(f)
+    isodd(p) || throw(DomainError(x, "inverse for x^$p is not defined at $x"))
+    copysign(abs(x)^inv(p), x)
 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
+    x ≥ zero(x) || throw(DomainError(x, "inverse for x^$p is not defined at $x"))
+    x^inv(p)
 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"))
-    end
+    isinteger(inv(p)) || throw(DomainError(x, "inverse for x^$p is not defined at $x"))
+    x^inv(p)
 end
 
 function invpow1(b::Real, x::Real)
+    # b < 0 should never happen in actual use: this check is done in inverse(f)
     if b ≥ zero(b) && x ≥ zero(x)
         log(b, x)
     else
@@ -44,14 +37,17 @@ function invpow1(b::Real, x::Real)
 end
 
 function invlog1(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
+    # exception may happen here: check cannot be done in inverse(f) because of log(Real, Complex)
+    b > zero(b) && b != one(b) || throw(DomainError(x, "inverse for log($b, x) is not defined at $x"))
+    b^x
 end
 invlog1(b, x) = b^x
 
+function invlog2(b::Real, x::Real)
+    # exception may happen here: check cannot be done in inverse(f) because of log(Complex, Real)
+    x > zero(x) && x != one(x) || throw(DomainError(x, "inverse for log($b, x) is not defined at $x"))
+    x^inv(b)
+end
 invlog2(b, x) = x^inv(b)
 
 

src/inverse.jl

@@ -147,17 +147,24 @@ 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(^), <:Integer}) = isodd(f.x) ? Base.Fix2(invpow2, f.x) : throw(DomainError(f.x, "Cannot invert x^$(f.x)"))
 inverse(f::Base.Fix2{typeof(invpow2)}) = Base.Fix2(^, f.x)
+
+inverse(f::Base.Fix1{typeof(^), <:Real}) = f.x > zero(f.x) ? Base.Fix1(invpow1, f.x) : throw(DomainError(f.x, "Cannot invert"))
 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(log)}) = f.x == one(f.x) ? throw(DomainError(f.x, "Cannot invert")) : 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(log)}) = f.x == one(f.x) ? throw(DomainError(f.x, "Cannot invert")) : Base.Fix2(invlog2, f.x)
 inverse(f::Base.Fix2{typeof(invlog2)}) = 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)
+
 inverse(f::Base.Fix2{typeof(fldmod)}) = Base.Fix2(invfldmod, f.x)
 inverse(f::Base.Fix2{typeof(invfldmod)}) = Base.Fix2(fldmod, f.x)
 

test/test_inverse.jl

@@ -46,7 +46,7 @@ InverseFunctions.inverse(f::Bar) = Bar(inv(f.A))
     x = rand()
     for f in (
             foo, inv_foo, log, log2, log10, log1p, sqrt,
-            Base.Fix2(^, rand()), Base.Fix2(^, rand([-10:-1; 1:10])), Base.Fix1(^, rand()), Base.Fix1(log, rand()), Base.Fix1(log, 1/rand()), Base.Fix2(log, rand()),
+            Base.Fix2(^, 3*rand() - 0.5), Base.Fix2(^, rand(float.([-10:-1; 1:10]))), Base.Fix1(^, rand()), Base.Fix1(log, rand()), Base.Fix1(log, 1/rand()), Base.Fix2(log, rand()),
         )
         InverseFunctions.test_inverse(f, x)
     end
@@ -80,7 +80,9 @@ InverseFunctions.inverse(f::Bar) = Bar(inv(f.A))
     @test_throws DomainError inverse(Base.Fix2(^, 0.5))(-5)
     @test_throws DomainError inverse(Base.Fix2(^, 0.51))(complex(-5))
     InverseFunctions.test_inverse(Base.Fix2(^, 0.5), complex(-5))
-    @test_throws DomainError inverse(Base.Fix2(^, 2))(-5)
+    @test_throws DomainError inverse(Base.Fix2(^, 2))
+    @test_throws DomainError inverse(Base.Fix2(^, -4))
+    InverseFunctions.test_inverse(Base.Fix2(^, 2.0), 4)
     @test_throws DomainError inverse(Base.Fix1(^, 2))(-5)
     @test_throws DomainError inverse(Base.Fix1(^, -2))(3)
     @test_throws DomainError inverse(Base.Fix1(^, -2))(3)