Fix irrational raised to a non-literal negative integer power

[lvlte]

Fixes #61284.

^(x::AbstractIrrational, y::Integer) was hitting (non-explicitly) ^(x::Number, p::Integer) = power_by_squaring(x,p).
Use float(x)^p which has a specialized method ^(x::Float64, n::Integer).

[adienes]

I wonder if instead, power_by_squaring should be changed so that rather than throw_domerr_powbysq(x, p) when p < 0, instead we return power_by_squaring(inv(x), -p)

[araujoms]

I think the most elegant solution would be to define the method

function ^(x::Union{AbstractFloat, AbstractIrrational, Rational}, p::Integer)
    if p < 0
        return power_by_squaring(inv(x), -p)
    else
        return power_by_squaring(x, p)
    end
end

It already exists for Rational, I don't know why it was omitted for AbstractFloat and AbstractIrrational.

[lvlte]

Yes, it would make the result "more consistent" with the literal power implementation :

julia/base/intfuncs.jl

Lines 479 to 489 in 984ad24

	@inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p}
	if p < 0
	if x isa BitInteger64
	f(Float64(x), p) # inv would cause rounding, while Float64^Integer is able to compensate the inverse
	else
	f(inv(x), -p)
	end
	else
	f(x, p)
	end
	end

And the suggestion by @adienes would also fix the issue for integers (I mean 2^-2 yields 0.25 but p=-2; 2^p throws a DomainError as for irrationals).

But I also noticed that the result when using the inverse - literal power and power_by_squaring - could be a bit off compared to float(x)^p, not sure if this always the case though (and if the contrary could be true) :

julia> x, p = pi, -10
(π, -10)

julia> Base.power_by_squaring(inv(x), -p)
1.0678279226861544e-5

julia> x^-10
1.067827922686154e-5

julia> float(x)^p
1.0678279226861537e-5

julia> Float64(big(x)^p)
1.0678279226861534e-5

This makes me think power_by_squaring was meant for integer base, which would explain why it throws when p is negative instead of using the inverse, and for rational the implementation mentioned by @araujoms does use the inverse probably because it is representable exactly.

[araujoms]

It's very much intentional that it doesn't work for integers; it would be otherwise type unstable. That it does work for literal powers is IMHO a bad decision, makes the language hacky and inconsistent. power_by_squaring is not necessarily meant for integers. It is used with more general objects like matrices, for example.

Maybe there is a good reason it wasn't done for AbstractIrrational, or maybe it was just an oversight.

[lvlte]

@araujoms This makes sense. Also maybe the reason AbstractFloat and AbstractIrrational not using power_by_squaring(inv(x), -p) as it's done with Rational is that the accuracy loss induced by the inversion makes the result worse than computing the arguments after promotion to floats (as in the example above).

[adienes]

actually, given the precedent

here

julia/base/irrationals.jl

Line 47 in 984ad24

promote_rule(::Type{<:AbstractIrrational}, ::Type{T}) where {T<:Real} = promote_type(Float64, T)

and here

julia/base/irrationals.jl

Line 197 in 984ad24

-(x::AbstractIrrational) = -Float64(x)

and here

julia/base/irrationals.jl

Line 59 in 984ad24

AbstractFloat(x::AbstractIrrational) = Float64(x)::Float64

maybe we just need ^(x::AbstractIrrational, p::Integer) = Float64(x)^p. or I guess for slightly more flexibility, float(x)^p. if that approach has less error than inv(x)^(-p)

[araujoms]

Then it's better to leave the PR as it is; promote is already taking it to Float64, I don't think it's necessary to hardcode Float64 again.

[araujoms]

@lvlte Ah, that's the reason! Indeed if we look at the implementation for Float64 we see that it computes the inverse, but also an error compensation term:

julia/base/special/pow.jl

Lines 122 to 128 in 984ad24

	if n < 0
	rx = inv(x)
	n==-2 && return rx*rx #keep compatibility with literal_pow
	isfinite(x) && (xnlo = -fma(x, rx, -1.) * rx)
	x = rx
	n = -n
	end

So my method was indeed a bad idea, it's better to stick to what you had already written.

[adienes]

I don't think just promote is as good though, since ^ promotes differently than other arithmetic operations. e.g. consider

julia> 0x02 * 3
6

julia> 0x02^3
0x08
```

[araujoms]

Argh, of course. And even worse, ^(x::AbstractIrrational, y::Integer) = ^(promote(x,y)...) will take it to ^(::Float64, Float64), which is not the method we want. I think then the best thing to do is ^(x::AbstractIrrational, y::Integer) = float(x)^y as you said.

I have to apologise for all this noise, I was trying to review the PR and take care of a baby simultaneously, and the result is that I performed poorly at both tasks.

[adienes]

odds of this breaking anything seem pretty low but I'd like to run nanosoldier just as a precaution

@nanosoldier runtests()