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

[stevengj]

This fixes #3024 by taking advantage of #20530 to transform x^-n to inv(x)^n for literal integer exponents n.

It means that things like 2^-2 now work (giving 0.25) in a type-stable way, you can use A^-1 for the inverse of a Matrix{Int}, etcetera.

(This pulls out the negative-exponent portions of PR #20637.)

[stevengj]

Note that this changes (and apparently greatly improves?) the code that is generated for e.g. f(x) = x^-2 when x::Float64. In Julia 0.6, I got:

julia> @code_llvm f(2.3)

define double @julia_f_60740(double) #0 !dbg !5 {
top:
  %1 = call double @llvm.pow.f64(double %0, double -2.000000e+00)
  %2 = fadd double %0, -2.000000e+00
  %notlhs = fcmp ord double %1, 0.000000e+00
  %notrhs = fcmp uno double %2, 0.000000e+00
  %3 = or i1 %notrhs, %notlhs
  br i1 %3, label %L12, label %if

if:                                               ; preds = %top
  call void @jl_throw(i8** inttoptr (i64 4624304072 to i8**))
  unreachable

L12:                                              ; preds = %top
  ret double %1
}

whereas now I get:

julia> @code_llvm f(2.3)
  %1 = fdiv double 1.000000e+00, %0
  %2 = fmul double %1, %1
  ret double %2
}

which looks a lot nicer. BenchmarkTools seems to be broken on 0.7 at the moment, but I can simulate the new behavior in 0.6 by g(x) = let y=inv(x); y*y; end, in which case I get:

julia> @btime f(2.3);
  37.351 ns (0 allocations: 0 bytes)

julia> @btime g(2.3);
  2.048 ns (0 allocations: 0 bytes)

which seems almost too good to be true. (They give exactly the same result f(2.3) == g(2.3).)

I really don't understand what the branch and jl_throw code are about in the f(x) code for Julia 0.6, though.

[KristofferC]

The branch is because the exponent 2 gets converted to 2.0 and then the function

@inline function ^(x::Float64, y::Float64)
    z = ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, y)
    if isnan(z) & !isnan(x+y)
        throw_exp_domainerror(x)
    end
    z
end

is called.

[stevengj]

Ah, thanks @KristofferC. Separate from this PR, maybe we should add a

@inline ^(x::Float64, y::Int) = ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, y)

method that skips the domain-error check for integer exponents (where I think it can never throw)?

(But I'm still confused: in what case does x^y throw a domain error for floating-point x and y? Oh, right, for negative values to fractional exponents.)

[KristofferC]

julia> -2.0^-2.3
-0.2030630990890589

julia> f(x, y) = x^y
f (generic function with 1 method)

julia> f(-2.0, -2.3)
ERROR: DomainError:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
 [1] nan_dom_err at ./math.jl:300 [inlined]
 [2] ^ at ./math.jl:699 [inlined]
 [3] f(::Float64, ::Float64) at ./REPL[2]:1

At http://en.cppreference.com/w/c/numeric/math/pow they say:

• If base is finite and negative and exponent is finite and non-integer, a domain error occurs and a range error may occur.
• If base is zero and exponent is zero, a domain error may occur.
• If base is zero and exponent is negative, a domain error or a pole error may occur.

[stevengj]

@KristofferC, we don't throw a domain error in the other two cases you mention. 0.0^0.0 gives 1.0 and 0.0^-1.0 gives Inf.

[KristofferC]

Well, it says "may occur"... :P

[stevengj]

@nanosoldier runbenchmarks(ALL, vs=":master")

[ararslan]

Nanosoldier doesn't work at the moment due to an issue with JLD on 0.7. I'm working on removing the dependency on JLD entirely, but until that's done, it's unlikely that any Nanosoldier run will succeed.

[nanosoldier]

Something went wrong when running your job:

NanosoldierError: failed to run benchmarks against primary commit: stored type BenchmarkTools.ParametersPreV006 does not match currently loaded type

Logs and partial data can be found here
cc @ararslan

[vtjnash]

Oh, I forgot to comment on this earlier: We weren't doing this in the past because of accuracy concerns (#8939 and #19872). We could easily have been doing this in LLVM and inference as an invisible optimization (instead of violating referential transparency) if that wasn't a concern.

[vtjnash]

this doesn't need to be @generated

[stevengj]

whoops, that was copy-and-paste typo, sorry.

[stevengj]

@vtjnash, we could only do this in LLVM for e.g. Float64, no? LLVM optimization wouldn't help with non-hardware types, and in particular it wouldn't allow 2^-3 etcetera.

As far as I can tell, #19872 is not an issue for the new literal_pow code.

Regarding #8939, you're right that there seems to be a slight degradation in accuracy for small powers, typically of 1ulp, but occasionally 2ulp for x^-2 and 4ulp for x^-5. So maybe we should revert to the old behavior for hardware floating-point types. (But use the inv transformation in @fastmath code?) Or maybe it is worth it to have a slight accuracy loss for the (really substantial?) performance gains.

[vtjnash]

It's also providing worse results with worse performance for non-hardware types like BigFloat, since we no longer fuse the operations, due to the catch-all behavior of this PR.

[stevengj]

With most types, e.g. complex numbers, the behavior is identical (since most types implement negative integer powers by inv(x)^n) or eliminates a DomainError (in cases that weren't supported before).

With BigFloat there is a very slight (3%) performance degradation because we don't re-use the temporary for inv(x), but to get good BigFloat performance we need to fix that in a more global way anyway. But sure, we should override the fallback for BigFloat.

Another way of saying it is that people implementing their own AbstractFloat type may want to override literal_pow to get optimal performance and accuracy. This is annoying. But there are very few such people, whereas there are a lot of people who are frustrated by domain errors for negative literal powers.

[StefanKarpinski]

I generally agree with @stevengj here, but I think @vtjnash's point is that we can just define ^ to be type-unstable and have return types like Int^Int :: Union{Int,Float64}. This can be made efficient for cases where the power is a literal by constant propagation. However, unless I'm missing something, this would introduce downstream type instability into any method calling Int^Int with non-literal exponent, which seems like an unacceptable cost. Sure, we could introduce more and cleverer optimizations to compensate for this problem, but we don't have those optimizations now and not introducing the performance trap in the first place seems like the more Julian way.

[fredrikekre]

xref #21014 (comment)