nthroot function

[stevengj]

@oscardssmith mentioned on discourse that IEEE 2018 recommends an nthroot(x, n) that computes the real nth-root of a real x, similar to the Matlab function nthroot and generalizing cbrt.

A possible implementation could be as simple as:

nthroot(x::Real, n::Integer) =
    isodd(n) || x ≥ 0 ? copysign(abs(x)^(1//n), x) :
    throw(DomainError("Exponentiation yielding a complex result requires a complex argument.  Replace nthroot(x, n) with complex(x)^(1//n)."))

though it's possible that we could have faster/more-accurate implementations in some special cases.

As discussed in #36534, it probably only makes sense to define this for x::Real arguments.

[mikmoore]

It seems that hermitian matrices should also have a compatible nthroot definition. That extension could be deferred to a follow-up PR, however.

EDIT: The discussion on discourse makes me think that we might actually want to provide a more general version of nthroot that extends to rational powers (proposed as realpow there, although that name wouldn't translate to hermitian matrices). nthroot(x,d)^n for x^(n//d) is probably passable for Real but would be rather inefficient for matrices.

[stevengj]

I would think that there should be an nthroot for any real matrix that returns a real result, similar to the discussion in JuliaLang/LinearAlgebra.jl#966? (Also an nthroot for Hermitian matrices that returns a Hermitian result.)

In principle, I agree that there might be a place for a realpow(x, p//n) that computes nthroot(x, n)^p, but note that this is very different from nthroot(x^p, n), which could be a bit confusing.

Alternatively, we could just have a method nthroot(x, n, p=1) to compute nthroot(x, n)^p, with a default definition nthroot(x::Real, n::Integer, p::Integer) = nthroot(x, n)^p?

[mikmoore]

I suppose triangular matrices with real diagonals also have real eigenvalues. So hermitian and real-diagonal triangular matrices are on the list. What other structures give rise to real eigenvalues? I don't know that I'd support having general matrices support such a function, since that could require hoping that numerical evaluation supports a nthroot-type root to avoid needing to throw an error.

Can you elaborate on the distinction between nthroot(x^p,n) vs nthroot(x,n)^p for reduced fractions p//n? I'm failing to see in what situations a difference arises (in infinite precision). If there is such a hazard, I'd be less supportive of rational powers.

[oscardssmith]

I think the problem is more that reducing the fraction changes the answer. E.G.

julia> ((complex(-1+0im))^(9))^(1/6)
0.8660254037844387 + 0.5im

julia> ((complex(-1+0im))^(3))^(1/2)
0.0 + 1.0im

[mikmoore]

Those two examples both have negative bases and even-denominator powers, so would result in DomainErrors with nthroot for any odd power-numerator (including 1) or any even power-numerator that failed to reduce the power-denominator to an odd. There's no real-valued solution for those examples regardless of the scaling of the exponent numerator/denominator. We've already established that nthroot(x,d) is distinct from ^(x,1//d) in returning real roots (or errors) rather than principle roots. So while nthroot(-1.0,6)^4 would error, any of nthroot((-1.0)^4,6), nthroot((-1.0)^2,3), or nthroot(-1.0,3)^2 would succeed and return an equivalent answer.

[Snimm]

Can I work on this?

[stevengj]

@Snimm, I think we need to first decide whether we want nthroot(x, p) or realpow(x, n//p)?

It sounds like people are leaning towards realpow(x, n//p)? Or should we have both?

• Advantage of nthroot: recommended by IEEE(?) and hence somewhat standardized(?), available in other languages (Matlab and R).
• Advantage of realpow: potentially slightly more efficient (especially if we extend to matrices, although even cbrt doesn't support matrices yet).
• Disadvantage of realpow: potentially confusing that realpow(-1.0, 4//6) is different than realpow(-1.0, 1//6)^4 (since the latter would error). On the other hand, in Julia 4//6 === 2//3 and realpow(-1.0, 2//3) == realpow(-1.0, 1//3)^2 == cbrt(-1.0)^2.

Or maybe we should have both? Or just start with nthroot?

[oscardssmith]

I'd start with nthroot

[matrixbot123]

I'd start with nthroot

what about making the cbrt function . if we use nthroot then cbrt will be less efficient . And the whole purpose of this issue was for writing cbrt function (kinda)

[oscardssmith]

we already have cbrt.

[matrixbot123]

we already have cbrt.

my mistake cbrt(A::AbstractMatrix{<:Real})

[StefanKarpinski]

Ultimately I think that having both makes sense. But implementation can be tackled starting with cbrt for matrices and nthroot. More general realpow can be done later seems like.

[mikmoore]

The established-ness of nthroot is certainly compelling. However, my reluctance to incorporate nthroot is that it is entirely (?) superseded by realpow. Once we add nthroot, we won't be able to remove it. We could remove it in v2.0, but we never would in practice.

Unless someone develops something profoundly more clever than the proposed implementation, the implementations are essentially identical so this feels rather wasteful. Meanwhile, there are some quantities that cannot be computed using just nthroot without carefully deciding which order to use among nthroot(x^m,n) vs nthroot(x,n)^m:

julia> (nextfloat(1.0)^(1//3))^10 # catastrophic loss of precision
1.0
julia> nextfloat(1.0)^(10//3)
1.0000000000000007
julia> (1e50^10)^(1//3) # overflow
Inf
julia> 1e50^(10//3)
4.641588833612859e166

Even once a viable order is decided, it still won't be faster or more accurate than the single exponentiation used in realpow.

Julia made the deliberate choice not to include the atan2 function (which is also recommended by IEEE754) because its functionality could be entirely captured under a 2-argument method of atan. I see the realpow/nthroot question as analogous.

P.S.
I don't have access to IEEE754, only to Wikipedia, but Wikipedia lists the operation as $x^\frac{1}{n}$ and does not specify a name. Does "nthroot" explicitly appear in the standard or is it merely a popular convention in languages and libraries? If the standard doesn't name it, then I'd lose my last reservations about abandoning nthroot.

[stevengj]

The only languages I can find with nthroot are Matlab and R. Mathematica calls it Surd. Fortran doesn't have it (see fortran-lang/stdlib#214).

So maybe we should just go with realpow, defined via e.g.:

realpow(x::Real, p::Rational) =
    isodd(p.den) || x ≥ 0 ? copysign(abs(x)^p, x) :
    throw(DomainError("Exponentiation yielding a complex result requires a complex argument.  Replace realpow(x, p) with complex(x)^p."))

(There is the question of what to do with -0.0 for even denominators. The above is consistent with sqrt(-0.0), which gives -0.0 according to IEEE 754.)

[mikmoore]

An issue with the name realpow is that the operation is well-defined for complex-valued Hermitian matrices (since the eigenvalues are real). nthroot makes no presumption about the real-ness of the inputs or outputs, but would be a poor name if rational powers are in-fact supported rather than just roots.

But realpow is an excellent name for use with real scalars and real triangular/diagonal matrices. So much so that I might accept the unfortunate confusion with Hermitian arguments (it's still correct in the eigenvalue sense). But if somebody can concoct a better name that reconciles this with Hermitian, I'd be all ears.

[stevengj]

The point is that realpow indicates which root is taken, whereas nthroot does not.

The matrix case is tricky (#47513) and I'm not sure that should be our main concern here. (For example, any real matrix has a real nth root when n is odd, which in my mind is a much more interesting use-case than applying this to complex-Hermitian matrices. The algorithms in the non-Hermitian matrix case are much harder, however, just as for matrix exponentials and square roots.)

[StefanKarpinski]

What about the name nthroot but with a three argument method that computes a rational power?

[stevengj]

Sure, why not.

[mikmoore]

I imagine the majority use of this function will, in practice, be to compute roots rather than general powers. So with that in mind I could live with nthroot and yet that function could support the bonus argument for a numerator. Although it does make the name a touch weird. I did appreciate the point that realpow did have the benefit of specifying the branch cut in the name.

[LilithHafner]

Would it be nthroot(base, denominator)/nthroot(base, numerator, denominator) or nthroot(base, denominator, numerator=1)?

If feels unfortunate to use one or two integers and argument order conventions to represent a rational rather than passing in an integer or rational directly.

[antoine-levitt]

Fwiw, I have a lot of trouble forcing my brain to parse nthroot as something else than n throot, and then I wonder what a throot is. I imagine it's even more confusing for people who are not very familiar with English. Realpow is a much better name.

[stevengj]

The big problem I have with nthroot is that, from the name, I expect it to be the same as x^(1//n) and work for any x. But we've learned to live with that for cbrt, so I can live with it either way.

[theWiseAman]

I can work on this one. Can anyone assign this to me? I have good idea of what to do here.
Plus how can I compile my function that I have written on local. I have already ran 'make' command for the entire julia directory but I don't want to do it for just one small change. Moreover how can I speed up makefile compilation. It's taking more than an hour 😥

[theWiseAman]

I wrote function for nthroot but this error is showing up

Am I missing something while compiling the makefile? I need a little help over the compilation of makefile.

[giordano]

Where did you write that? If it's inside the Base module then you'll also need to export it to call the function unqualified, otherwise you can always call it by prefixing the module it lives in, Base.nthroot or wherever that is. Also, for quick tests you don't need to define the function inside Julia source code and recompile it every time, that'd be a pretty slow development cycle, you can do that directly in the REPL

[theWiseAman]

Where did you write that? If it's inside the Base module then you'll also need to export it to call the function unqualified, otherwise you can always call it by prefixing the module it lives in, Base.nthroot or wherever that is. Also, for quick tests you don't need to define the function inside Julia source code and recompile it every time, that'd be a pretty slow development cycle, you can do that directly in the REPL

@giordano Ya its in base/complex.jl file. Ok let me see how to export it. Base.nthroot is working.
So, what should be the efficient development procedure. For example, if I modified the code a little bit and want to see the results in the julia shell, so what should be my immediate steps. Should I just run make clean then make?
I would like to document the answer of this for future newbies like me.

[oscardssmith]

use Revise. you can just Revise.track(Base) and then just saving the file you are editing will update the Julia repl

[theWiseAman]

use Revise. you can just Revise.track(Base) and then just saving the file you are editing will update the Julia repl

@oscardssmith Ya I need to read docs properly 😥. I will see how to do that. Can you check the Draft PR and give suggestions? Also is my function exported correctly?

PS: Thanks a bunch mate :)
PS: For future reference
Julia REPL means the julia shell executable. For the first time we need to install Revise package by following this Easy Revise Installation. Then use this guide Revise Guide.

[LilithHafner]

I can work on this one. Can anyone assign this to me?

Thanks! There are some tricky unresolved design decisions here which we will need to resolve before merging a PR, so this might not be the best first issue. Nevertheless, you are certainly welcome to make a PR! You don't have to be assigned to an issue to contribute a PR addressing it.

Can you check the draft PR and give suggestions? Also is my function exported correctly?

It would be easier to find if you linked it with ...check the [draft PR](47860) and.... It seems like you've exported it properly but I'm not in a position to comment on whether the semantics are right.

[theWiseAman]

so this might not be the best first issue

Don't worry I got this now

Thanks! There are some tricky unresolved design decisions here which we will need to resolve before merging a PR

Can you comment more specifically with regards to my Draft PR (#47860)? What changes are needed in the function I wrote?

[11happy]

@stevengj I am interested to work on this,but based on the discussions and work above I am quite confused what to implement realpow or nthroot , for Real + Complex numbers or only for Real, Can you please guide me a bit.
Thank you

[stevengj]

See my comment here: #47860 (comment)

I think it's pretty clear that we want something only for Real, and only returns the real root following the IEEE functionality. It's less clear whether we want nthroot(x::Real, n::Integer), realpow(x::Real, p::Rational), or nthroot(x::Real, n::Integer, p::Integer), or …