Rational power of BigFloat should use specialised MPFR function

[dpsanders]

big(1.3)^(1//7) currently converts the rational to BigFloat before computing the power, thereby losing accuracy.

It should instead use the specialised MPFR function mpfr_rootn_si.

cc @oscardssmith

[dpsanders]

Apparently the code is still in the library

[dpsanders]

julia> x = interval(8.673020346900622e8, 8.673020346900623e8)

julia> IntervalArithmetic._rootn_round(sup(x), 8, RoundUp)
13.100000000000001

[dpsanders]

At some point we had an nthroot function. This should maybe be renamed rootn and reinstated?

[dpsanders]

(The rootn function is required recommended by the 1788 standard in section 10.6.)

[dpsanders]

My bad, I see we have the rootn function already... In that case we should just make rational power use that function!

[dpsanders]

That works if the numerator == 1, so the rational is 1 // n.
If the rational is x ^ (p // q) then we could do nthroot(x, q)^p. This will not be guaranteed to be correctly rounded.

[oscardssmith]

Note that nthroot(x, q) is going to lose a lot of accuracy for x close to 0.

[ararslan]

Would something like this be a good first-pass solution?

diff --git a/base/mpfr.jl b/base/mpfr.jl
index ed3ea5937c..7833b65284 100644
--- a/base/mpfr.jl
+++ b/base/mpfr.jl
@@ -733,6 +733,25 @@ function ^(x::BigFloat, y::BigInt)
     return z
 end

+_rootn(z, x, n::ClongMax)  = ccall((:mpfr_rootn_si, libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Clong,  MPFRRoundingMode), z, x, n, rounding_raw(BigFloat))
+_rootn(z, x, n::CulongMax) = ccall((:mpfr_rootn_ui, libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Culong, MPFRRoundingMode), z, x, n, rounding_raw(BigFloat))
+
+function ^(x::BigFloat, y::Rational{<:Union{ClongMax,CulongMax}})
+    # If `y` is of the form `±1//n`, we can use MPFR's builtin nth root without needing to
+    # convert the `Rational` to a `BigFloat`. The `copysign` is because `mpfr_rootn_*` takes
+    # a negative power to mean `-1//n` but `Rational{<:Signed}` uses the numerator to track
+    # the sign.
+    if isone(abs(numerator(y)))
+        z = BigFloat()
+        _rootn(z, x, copysign(denominator(y), numerator(y)))
+        return z
+    else
+        # TODO: Determine whether there are cases where we should prefer doing something
+        # like `x^(p//q) = rootn(x, q)^p`, which will have accuracy issues for `x` near 0
+        return x^BigFloat(y)
+    end
+end
+
 ^(x::BigFloat, y::Integer)  = typemin(Clong)  <= y <= typemax(Clong)  ? x^Clong(y)  : x^BigInt(y)
 ^(x::BigFloat, y::Unsigned) = typemin(Culong) <= y <= typemax(Culong) ? x^Culong(y) : x^BigInt(y)

[oscardssmith]

Actually thinking about this more, as long as x^p is neither 0 or inf, nthroot(x^p, q) will be pretty accurate.