make x^-n equivalent to inv(x)^n for literal n

NEWS.md

@@ -259,6 +259,10 @@ Library improvements
     If this argument is used they return a string consisting of first/last `nchar`
     characters from the original string ([#23960]).
 
+  * Expressions `x^-n` where `n` is an *integer literal* now correspond to `inv(x)^n`.
+    For example, `x^-1` is now essentially a synonym for `inv(x)`, and works
+    in a type-stable way even if `typeof(x) != typeof(inv(x))` ([#24240]).
+
   * The functions `nextind` and `prevind` now accept `nchar` argument that indicates
     the number of characters to move ([#23805]).
 

base/fastmath.jl

@@ -93,6 +93,9 @@ const rewrite_op =
 function make_fastmath(expr::Expr)
     if expr.head === :quote
         return expr
+    elseif expr.head == :call && expr.args[1] == :^ && expr.args[3] isa Integer
+        # mimic Julia's literal_pow lowering of literal integer powers
+        return Expr(:call, :(Base.FastMath.pow_fast), make_fastmath(expr.args[2]), Val{expr.args[3]}())
     end
     op = get(rewrite_op, expr.head, :nothing)
     if op !== :nothing
@@ -263,6 +266,8 @@ end
 
 pow_fast(x::Float32, y::Integer) = ccall("llvm.powi.f32", llvmcall, Float32, (Float32, Int32), x, y)
 pow_fast(x::Float64, y::Integer) = ccall("llvm.powi.f64", llvmcall, Float64, (Float64, Int32), x, y)
+pow_fast(x::FloatTypes, ::Val{p}) where {p} = pow_fast(x, p) # inlines already via llvm.powi
+@inline pow_fast(x, v::Val) = Base.literal_pow(^, x, v)
 
 sqrt_fast(x::FloatTypes) = sqrt_llvm(x)
 

base/intfuncs.jl

@@ -234,6 +234,18 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}}
 @inline literal_pow(::typeof(^), x::HWNumber, ::Val{2}) = x*x
 @inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x
 
+@inline @generated function literal_pow(f::typeof(^), x, ::Val{p}) where {p}
+    if p < 0
+        :(literal_pow(^, inv(x), $(Val{-p}())))
+    else
+        :(f(x,$p))
+    end
+end
+
+# note: it is tempting to add optimized literal_pow(::typeof(^), x, ::Val{n})
+#       methods here for various n, but this easily leads to method ambiguities
+#       if anyone has defined literal_pow(::typeof(^), x::T, ::Val).
+
 # b^p mod m
 
 """

base/linalg/generic.jl

@@ -825,6 +825,11 @@ function pinv(v::AbstractVector{T}, tol::Real=real(zero(T))) where T
     return res
 end
 
+# this method is just an optimization: literal negative powers of A are
+# already turned by literal_pow into powers of inv(A), but for A^-1 this
+# would turn into inv(A)^1 = copy(inv(A)), which makes an extra copy.
+@inline Base.literal_pow(::typeof(^), A::AbstractMatrix, ::Val{-1}) = inv(A)
+
 """
     \\(A, B)
 

base/mathconstants.jl

@@ -85,6 +85,7 @@ catalan
 for T in (Irrational, Rational, Integer, Number)
     Base.:^(::Irrational{:ℯ}, x::T) = exp(x)
 end
+@generated Base.literal_pow(::typeof(^), ::Irrational{:ℯ}, ::Val{p}) where {p} = exp(p)
 
 Base.log(::Irrational{:ℯ}) = 1 # use 1 to correctly promote expressions like log(x)/log(ℯ)
 Base.log(::Irrational{:ℯ}, x::Number) = log(x)

doc/src/manual/faq.md

@@ -274,13 +274,6 @@ ERROR: DomainError with -2.0:
 sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
 Stacktrace:
 [...]
-
-julia> 2^-5
-ERROR: DomainError with -5:
-Cannot raise an integer x to a negative power -5.
-Make x a float by adding a zero decimal (e.g., 2.0^-5 instead of 2^-5), or write 1/x^5, float(x)^-5, or (x//1)^-5
-Stacktrace:
-[...]
 ```
 
 This behavior is an inconvenient consequence of the requirement for type-stability.  In the case
@@ -296,9 +289,6 @@ your willingness to accept an *output type* in which the result can be represent
 ```jldoctest
 julia> sqrt(-2.0+0im)
 0.0 + 1.4142135623730951im
-
-julia> 2.0^-5
-0.03125
 ```
 
 ### Why does Julia use native machine integer arithmetic?

test/complex.jl

@@ -784,8 +784,10 @@ end
 
 @test complex(1//2,1//3)^2 === complex(5//36, 1//3)
 @test complex(2,2)^2 === complex(0,8)
-@test_throws DomainError complex(2,2)^(-2)
-@test complex(2.0,2.0)^(-2) === complex(0.0, -0.125)
+let p = -2
+    @test_throws DomainError complex(2,2)^p
+end
+@test complex(2,2)^(-2) === complex(2.0,2.0)^(-2) === complex(0.0, -0.125)
 
 @test complex.(1.0, [1.0, 1.0]) == [complex(1.0, 1.0), complex(1.0, 1.0)]
 @test complex.([1.0, 1.0], 1.0) == [complex(1.0, 1.0), complex(1.0, 1.0)]

test/fastmath.jl

@@ -208,3 +208,7 @@ end
     @fastmath a[idx...] += b[idx...]
     @test a == b
 end
+
+@testset "literal powers" begin
+    @test @fastmath(2^-2) == @fastmath(2.0^-2) == 0.25
+end

test/linalg/dense.jl

@@ -701,7 +701,9 @@ end
     @testset "Tests for $elty" for elty in (Int128, Int16, Int32, Int64, Int8,
                                             UInt128, UInt16, UInt32, UInt64, UInt8,
                                             BigInt)
-        @test_throws DomainError elty[1 1;1 0]^-2
+        info("Testing $elty")
+        @test elty[1 1;1 0]^-1 == [0  1;  1 -1]
+        @test elty[1 1;1 0]^-2 == [1 -1; -1  2]
         @test (@inferred elty[1 1;1 0]^2) == elty[2 1;1 1]
         I_ = elty[1 0;0 1]
         @test I_^-1 == I_

test/linalg/symmetric.jl

@@ -262,18 +262,10 @@ end
                 @testset "pow" begin
                     # Integer power
                     @test (asym)^2   ≈ (Symmetric(asym)^2)::Symmetric
-                    if eltya <: Integer && !isone(asym) && !isone(-asym)
-                        @test_throws DomainError (asym)^-2
-                    else
-                        @test (asym)^-2  ≈ (Symmetric(asym)^-2)::Symmetric
-                    end
+                    @test (asym)^-2  ≈ (Symmetric(asym)^-2)::Symmetric
                     @test (aposs)^2  ≈ (Symmetric(aposs)^2)::Symmetric
                     @test (aherm)^2  ≈ (Hermitian(aherm)^2)::Hermitian
-                    if eltya <: Integer && !isone(aherm) && !isone(-aherm)
-                        @test_throws DomainError (aherm)^-2
-                    else
-                        @test (aherm)^-2 ≈ (Hermitian(aherm)^-2)::Hermitian
-                    end
+                    @test (aherm)^-2 ≈ (Hermitian(aherm)^-2)::Hermitian
                     @test (apos)^2   ≈ (Hermitian(apos)^2)::Hermitian
                     # integer floating point power
                     @test (asym)^2.0   ≈ (Symmetric(asym)^2.0)::Symmetric

test/math.jl

@@ -619,14 +619,18 @@ end
     end
 end
 
-@testset "issues #3024, #12822" begin
-    @test_throws DomainError 2 ^ -2
+@testset "issues #3024, #12822, #24240" begin
+    p2 = -2
+    p3 = -3
+    @test_throws DomainError 2 ^ p2
+    @test 2 ^ -2 == 0.25 == (2^-1)^2
     @test_throws DomainError (-2)^(2.2)
     @test_throws DomainError (-2.0)^(2.2)
-    @test_throws DomainError false ^ -2
-    @test 1 ^ -2 === (-1) ^ -2 === 1
-    @test (-1) ^ -3 === -1
-    @test true ^ -2 === true
+    @test_throws DomainError false ^ p2
+    @test false ^ -2 == Inf
+    @test 1 ^ -2 === (-1) ^ -2 == 1 ^ p2 === (-1) ^ p2 === 1
+    @test (-1) ^ -1 === (-1) ^ -3 == (-1) ^ p3 === -1
+    @test true ^ -2 == true ^ p2 === true
 end
 
 @testset "issue #13748" begin

test/numbers.jl

@@ -2944,16 +2944,19 @@ Base.literal_pow(::typeof(^), ::PR20530, ::Val{p}) where {p} = 2
     @test [x,x,x].^2 == [2,2,2]
     for T in (Float16, Float32, Float64, BigFloat, Int8, Int, BigInt, Complex{Int}, Complex{Float64})
         for p in -4:4
-            if p < 0 && real(T) <: Integer
-                @test_throws DomainError eval(:($T(2)^$p))
-            else
-                v = eval(:($T(2)^$p))
-                @test 2.0^p == T(2)^p == v
+            v = eval(:($T(2)^$p))
+            @test 2.0^p == v
+            if p >= 0 || T == float(T)
+                @test v == T(2)^p
                 @test v isa T
+            else
+                @test v isa float(T)
             end
         end
     end
     @test PR20889(2)^3 == 5
+    @test [2,4,8].^-2 == [0.25, 0.0625, 0.015625]
+    @test ℯ^-2 == exp(-2) ≈ inv(ℯ^2) ≈ (ℯ^-1)^2 ≈ sqrt(ℯ^-4)
 end
 module M20889 # do we get the expected behavior without importing Base.^?
     using Test

test/replutil.jl

@@ -256,7 +256,7 @@ struct TypeWithIntParam{T <: Integer} end
 let undefvar
     err_str = @except_strbt sqrt(-1) DomainError
     @test contains(err_str, "Try sqrt(Complex(x)).")
-    err_str = @except_strbt 2^(-1) DomainError
+    err_str = @except_strbt 2^(1-2) DomainError
     @test contains(err_str, "Cannot raise an integer x to a negative power -1")
     err_str = @except_strbt (-1)^0.25 DomainError
     @test contains(err_str, "Exponentiation yielding a complex result requires a complex argument")