improvements to accuracy/performance for float^integer

[stevengj]

This PR makes a few improvements to accuracy and/or performance for float^integer exponentiation.

• Float64^integer and Float32^integer no longer check for DomainErrors, since these can only occur for fractional powers (make x^-n equivalent to inv(x)^n for literal n #24240 (comment)). This seems to greatly improve performance in some cases (by a factor of 20+) on my machine, and the performance advantage of inv(x)^n over x^-n seems to be mostly gone.

• float^-n no longer uses the inv(x)^n fallback for literal n, since that degrades accuracy slightly (make x^-n equivalent to inv(x)^n for literal n #24240 (comment))

[stevengj]

CI failures seem unrelated.

[stevengj]

Ping @StefanKarpinski.

[StefanKarpinski]

Looks good to me but I'm not sure I'm the best person to review this. Maybe @vtjnash for the LLVM bit and @simonbyrne for the math bit?

[simonbyrne]

Is it worth adding @inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{-1}) = 1/x as well?

[vtjnash]

Why are these ambiguous? LLVM already knows about how to do this rewrite, when it can show that it won't cause accuracy issues. That at least seemed to be the conclusion in #19890 (otherwise, I would have also just altered the flag to let LLVM know that it was allowed to do this rewrite).

[stevengj]

Without these methods, literal_pow(::typeof(^), x::AbstractFloat, ::Val{p}) is ambiguous with literal_pow(::typeof(^), x::HWNumber, ::Val{2}) etc.

[stevengj]

@simonbyrne, that couldn't hurt, I guess — might as well add it.

[stevengj]

@vtjnash, actually, it seems like there is no longer an ambiguity — there used to be an ambiguity when I was defining literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{p}), but it looks like it may be fine now that it is AbstractFloat.

[vtjnash]

Hm, ok. It still seems unnecessary to be defining a special function overlay to be doing something that LLVM is doing already (e.g., it already rewrites pow(x, 2) as x*x whenever possible). The same is true of -1 method proposed above: that optimization is already implemented transparently in the backend.

[stevengj]

LLVM only does that conversion for types that it knows about. e.g. it won't do it for Complex{Float64} or DecimalFloat32.

[vtjnash]

It seems like you should have to opt-in to getting a less accurate answer that violates referential transparency, rather than opt-out of it.

[stevengj]

For most types, it's not less accurate. For most floating-point-based types (complex numbers, matrices, etc.), negative integer powers are already computed by inv(x)^-n. And for many other types (e.g. Int), negative literal powers wouldn't be defined at all without this transformation.

Can we not re-litigate the whole literal_pow debate in this PR? The question is whether this PR is an improvement on the current literal_pow rules.

[vtjnash]

Why are these ambiguous? LLVM already knows about how to do this rewrite, when it can show that it won't cause accuracy issues. That at least seemed to be the conclusion in #19890 (otherwise, I would have also just altered the flag to let LLVM know that it was allowed to do this rewrite).

[stevengj]

Bump.

[stevengj]

Any chance of merging this? It seems like a strict improvement over the current situation, and CI passes except for the usual AppVeyor hiccups.

[iblislin]

master build failed after this PR merged:

Error in testset numbers:
Test Failed at /home/julia/ci/worker/11rel-amd64/build/test/numbers.jl:2971
  Expression: 0.09496527f0 ^ -2 === 110.88438f0
   Evaluated: 110.884384f0 === 110.88438f0

seems only happened on amd64 build

Ref:

• https://ci.appveyor.com/project/JuliaLang/julia/build/1.0.21893/job/etytka7vk8p7kl2w
• https://julia.iblis.cnmc.tw/#/builders/1/builds/5028
• https://circleci.com/gh/JuliaLang/julia/10891

[fredrikekre]

Should have rerun CI here since it was over 2 weeks since the last time.

[StefanKarpinski]

Sorry, my bad.

[KristofferC]

So what do we do? Quick revert or live with it until someone fixes it?

[StefanKarpinski]

@stevengj – do you think you could fix this soon or should I revert this?

[StefanKarpinski]

Here's the revert PR: #24930

[simonbyrne]

It appears to be due to the use of LLVM intrinsics. I think we've had trouble with those before (e.g. #2741 & #8939).

[stevengj]

I don't really understand why there would be a problem with the LLVM intrinsics here. The ccall is virtually identical to the one in ^(x::Float64, y::Float64), just with the subsequent isnan check removed.

[stevengj]

I'm guessing that the code is perfectly fine, but that 0.09496527f0 ^ -2 === 110.88438f0 test is too stringent? It looks like on amd64 it is actually getting a more accurate result, which should be fine…

(Is it something about extended-precision mode on Linux?)