Fix Float32/16 raised to integer typemin

[kuszmaul]

Fixes #57464

The code for x^n where x::Float32 and n::Int previously failed for Float32(1.1)^typemin(Int) because it would reduce the problem to inv(x)^-n. That works fine unless n is typemin, in which case n==-n.

This PR makes a special case for n==typemin to effectively compute (x^(n/2))^2. Just in case n is also odd, we do x^cld(n,2) * x^fld(n,2).

[LilithHafner]

Thanks for working on this! You picked a great first issue to contribute based on.

[oscardssmith]

This looks like a functional fix, but I think it might be better to use the same approach that #53967 adds where we use a floating point algorithm for large powers (especially since it would turn this PR into a performance improvement rather than a regression. Specifically, I think the following would work well. (testing now)

@constprop :aggressive @inline function ^(x::Float32, n::Integer)
    n = clamp(n, Int64)
    n == 0 && return one(x)
    if use_power_by_squaring(n)
        n < 0 && return oftype(x, Base.power_by_squaring(inv(widen(x)), -n))
        return oftype(x, Base.power_by_squaring(inv(widen(x)), n))
    else
        s = ifelse(x < 0 && isodd(n), -1f0, 1f0)
        x = abs(x)
        y = float(n)
        return copysign(Float32(exp2(log2(Float64(x)*y))), s)
    end
end

[kuszmaul]

It does look like a performance improvement (once you get y outside of the log2. It should be exp2(log2(Float64(x))*y)
And it seems not to suffer from any errors in precision. You can just use pow_body.

[kuszmaul]

Although I don't think this PR is a regression. It runs the same code that it always ran for cases except typemin, which previously crashed.

[oscardssmith]

I meant performance regression. This code is pretty fast to the point that extra branches can have a significant speed impact. This way the branch is still there, but it is providing a speedup for lots of large powers as well.

[kuszmaul]

The branch is highly predictable (it's almost never taken), so it probably doesn't have any performance impact. This path already has many branches that aren't so predictable. Is there some evidence of a performance regression?

The x^y=exp(log(x)*y) approach does seem faster, however, and seems likely to be just as accurate with widening. Do you plan to prepare a PR for it?

[mikmoore]

Whatever we do needs to give (-1.0f0) ^ (2^60+1) == - (-1.0f0) ^ (2^60+2), so be careful with converting the exponent to a float.

This issue looks like it could be solved with the one-line change of replacing -n with Base.uabs(n). This is exactly what Base.uabs is for.

[oscardssmith]

That case is already handled (since s = ifelse(x < 0 && isodd(n), -1f0, 1f0 will flip the sign based on n and log2(x) is 0 so y won't factor into the equation at all.

[LilithHafner]

@oscardssmith, do you plan to prepare a PR with your suggestion, as @kuszmaul asked? Or would you rather put that implementation into this PR?

[oscardssmith]

oops, sorry for forgetting this. I'm happy to make the PR.

[KlausC]

Why not simply for all n < 0:

    i = inv(widen(x))
    return oftype(x, Base.power_by_squaring(i, -(n+1)) * i)

[KlausC]

cases n == 0, n == 1and n == -1 perform unnecessary conversions (widen and recast).

[KlausC]

The fix seems more complex than necessary. The cases abs(n) <= 2 should be special cased IMO.

[KlausC]

The same is true, if you use for example Int8 or Int16. So why define Int3 ?

[oscardssmith]

Is this PR not fixed by #57829?

[KlausC]

Is this PR not fixed by #57829?

Yes, it is superseded by that - I left similar comments there about the test cases. Unfortunately I saw this PR first.

[KlausC]

The #57829 is much better, of course, ant this PR should be closed IMO.