Resolves #47565 and #36534. nthroot() for complex variables implemented.

[theWiseAman]

So I have written the function nthroot. It returns all the roots in 1D array be it real or complex.

You can see in the above picture that not only it matches sqrt() but also improves it, as the sqrt() reports only one root instead of two. First 4 lines are the examples that the my code is working correctly.

Though I see two points of improvement.

1. First is related to the precision handling of float numbers (see for example it shows 10^(-16) rather than 0, also there is some slight fluctuations for 0.5).
2. Second, I want that nthroot(-1, 3) and nthroot(-1+0im, 3) to return same thing (will it be a good idea?). Or should I implement nthroot(x::Real, n::Int) so as to return a real value from the roots array if it exits, else throw an error like below
   

Feedback would be appreciated!! Also is my function exported correctly?

[bramtayl]

This seems potentially possible to do as an iterator to avoid the allocations?

[theWiseAman]

@bramtayl Is this what you asked for?

diff --git a/base/complex.jl b/base/complex.jl
index 27fc62dfff..67e9721126 100644
--- a/base/complex.jl
+++ b/base/complex.jl
@@ -561,7 +561,7 @@ end
 #     return Complex(abs(iz)/r/2, copysign(r,iz))
 # end
 
-function nthroot(z::Complex, n::Int)
+function nthroot(z::Complex, n::Integer)
     z = float(z)
     x, y = reim(z)
     if x==y==0
@@ -569,9 +569,9 @@ function nthroot(z::Complex, n::Int)
     end
     r = abs(z)
     θ = angle(z)
-    roots = Array{Complex}(undef, n)
+    roots = ComplexF64[]
     for k in 1:n
-        roots[k] = r^(1/n)*(cos((θ+2*pi*(k-1))/n) + sin((θ+2*pi*(k-1))/n)im)
+        push!(roots, ComplexF64(r^(1/n)*(cis((θ+2*pi*(k-1))/n))))
     end
     return roots
 end

Tries to solve issues #36534 and #47565

[bramtayl]

Not quite. I was thinking something more along the lines of RootArray <: AbstractArray. You would have to implement the abstract array interface here: https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array

[theWiseAman]

Sample results for nthroot()

Edit: Handling of complex domain errors after build fixes

[theWiseAman]

Though I see two points of improvement.

1. First is related to the precision handling of float numbers (see for example it shows 10^(-16) rather than 0, also there is some slight fluctuations for 0.5).

Improved my code along this idea.

2. Second, I want that nthroot(-1, 3) and nthroot(-1+0im, 3) to return same thing (will it be a good idea?). Or should I implement nthroot(x::Real, n::Int) so as to return a real value from the roots array if it exits, else throw an error like below

Overloaded nthroot(): nthroot(x::Real, n::Integer) and nthroot(z::Complex, n::Integer)

[theWiseAman]

Not quite. I was thinking something more along the lines of RootArray <: AbstractArray. You would have to implement the abstract array interface here: https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array

@bramtayl Give me some time to look into it. Rest is my function correct mathematically?
Just One thing remains, when x::Real, n::even and x > 0 my function should return two roots I suppose.

[LilithHafner]

I like returning a RootVector <: AbstractVector or perhaps accepting a third argument to nthroot that specifies which root so long as it is possible to make first(nthroots(x, n)) or nthroot(x, n, 1) run as fast as the less general nthroot(x, n) which only returns one root.

Even when returning all roots, there is a question about where to put the branch point. Where does this PR put it? i.e. for which x and n is the function first(nthroots(x, n)) discontinuous?

As an aside you can post code like this

```julia
julia> sqrt(2)
1.4142135623730951
```

julia> sqrt(2)
1.4142135623730951

which is better for folks with weak internet connections, renders the same way no matter who is posting it (unlike screenshots which depend on the editor settings of the person posting them) and more importantly, can be copied.

[theWiseAman]

I like returning a RootVector <: AbstractVector

I think now I am getting the rough idea of what is required here. Lemme refer the Docs

[bramtayl]

This example might be helpful:

julia> struct SquaresVector <: AbstractArray{Int, 1}
           count::Int
       end

julia> Base.size(S::SquaresVector) = (S.count,)

julia> Base.IndexStyle(::Type{<:SquaresVector}) = IndexLinear()

julia> Base.getindex(S::SquaresVector, i::Int) = i*i

[theWiseAman]

@LilithHafner @bramtayl
Is this what is required?

struct RootsVector{T<:Complex} <: AbstractVector{T}
    z::T
    n::Int
end

Base.size(S::RootsVector) = S.n
Base.IndexStyle(S::RootsVector) = IndexLinear()
Base.getindex(S::RootsVector, i::Int) = abs(S.z)^(1/S.n)*(Float16(cos((angle(S.z)+2*pi*(i-1))/S.n)) + Float16(sin((angle(S.z)+2*pi*(i-1))/S.n))im)

But I am unable to figure out two things:-

1. How to export RootsVector from base/complex.jl?
2. I guess I am missing something because of which Julia REPL is throwing print errors as below, though roots gets all the values filled in its indices. (Edit: Fixed)

julia> roots = RootsVector(1+0im, 3)
ERROR: UndefVarError: `RootsVector` not defined
Stacktrace:
 [1] top-level scope
   @ REPL[1]:1

julia> roots = Base.RootsVector(1+0im, 3)
1×1×1 Base.RootsVector{Complex{Int64}} with indices 1×2×3:
Error showing value of type Base.RootsVector{Complex{Int64}}:
ERROR: MethodError: no method matching unitrange(::Int64)

julia> for i in 1:3
           println(roots[i])
       end
1.0 + 0.0im
-0.5 + 0.8662109375im
-0.5 - 0.8662109375im

Sample Example that it is returning correct values

julia> roots = Base.RootsVector(1+0im, 4)
1×1×1×1 Base.RootsVector{Complex{Int64}} with indices 1×2×3×4:
Error showing value of type Base.RootsVector{Complex{Int64}}:
ERROR: MethodError: no method matching unitrange(::Int64)

julia> for i in 1:4
           println(roots[i])
       end
1.0 + 0.0im
0.0 + 1.0im
-1.0 + 0.0im
-0.0 - 1.0im

[LilithHafner]

Base.size(S::RootsVector) = S.n

size Returns a tuple containing the dimensions of A

julia> size([1,2,3])
(3,)

julia> size(RootsVector(1+0im, 3))
3

[theWiseAman]

@LilithHafner @bramtayl

But I am unable to figure out two things:-

1. How to export RootsVector from base/complex.jl?

How to export it though?

2. I guess I am missing something because of which Julia REPL is throwing print errors as below, though roots gets all the values filled in its indices.

This is fixed now.

julia> Base.RootsVector(1+0im, 3)
3-element Base.RootsVector{Complex{Int64}}:
  1.0 + 0.0im
 -0.5 + 0.8662109375im
 -0.5 - 0.8662109375im

julia> Base.RootsVector(1+0im, 4)
4-element Base.RootsVector{Complex{Int64}}:
  1.0 + 0.0im
  0.0 + 1.0im
 -1.0 + 0.0im
 -0.0 - 1.0im

Even when returning all roots, there is a question about where to put the branch point. Where does this PR put it? i.e. for which x and n is the function first(nthroots(x, n)) discontinuous?

I didn't quite understood what you mean by branch point (some examples would be preferable), but I am looking into it how to incorporate it. Shouldn't branch point be zero(1im) i.e. 0 + 0im for w = z^(1/n) where n::Integer?

julia> sqrt(0+0im)
0.0 + 0.0im

julia> Base.RootsVector(0+0im, 4)
4-element Base.RootsVector{Complex{Int64}}:
  0.0 + 0.0im
  0.0 + 0.0im
 -0.0 + 0.0im
 -0.0 - 0.0im

julia> Base.RootsVector(0+0im, 3)
3-element Base.RootsVector{Complex{Int64}}:
  0.0 + 0.0im
 -0.0 + 0.0im
 -0.0 - 0.0im

[LilithHafner]

How to export it though?

You can export it with export RootVector if the interface your prefer is RootVector(x, n).

some examples would be preferable

sqrt has a branch cut along the negative real line. See also https://en.wikipedia.org/wiki/Branch_point#Branch_cuts

julia> x = -1
-1

julia> Δ = 1e-9im
0.0 + 1.0e-9im

julia> sqrt(x + Δ) - sqrt(x - Δ)
0.0 + 2.0im

[mcabbott]

It seems quite surprising for nthroot(8, 3) to be a number and nthroot(8+0im, 3) a vector. If both are useful, maybe the one which returns n roots should be nthroots(x, n) plural?

For implementing this, using Float16 to round will give terrible accuracy except in easy cases (like those printed above). You may want some strategy of calculating in a higher precision than the final answer, like using Float64 for Float32 answers... and perhaps Base.TwicePrecision for Float64, but that's a whole can of worms.

[thchr]

Examples of problems where all roots are actually needed would be interesting to understand the use case for a multivalued nthroot.
They are all related by roots of unity, so on the face of it, the utility of having them all is quite small?

[theWiseAman]

For implementing this, using Float16 to round will give terrible accuracy except in easy cases (like those printed above). You may want some strategy of calculating in a higher precision than the final answer, like using Float64 for Float32 answers... and perhaps Base.TwicePrecision for Float64, but that's a whole can of worms.

It was representing zero like below that's why I used Float16, but I agree with your point. Is it fine that it shows zero like below?

julia> roots = RootsVector(1+0im, 4)
4-element RootsVector{Complex{Int64}}:
                     1.0 + 0.0im
   6.123234262925839e-17 + 1.0im
                    -1.0 + 1.2246468525851679e-16im
 -1.8369701465288538e-16 - 1.0im

[theWiseAman]

They are all related by roots of unity, so on the face of it, the utility of having them all is quite small?

Do we have a function which returns roots of unity? I am writting this to generalize for complex numbers beyond just unity. I think n roots might be useful, so we should have that option to extract it.

[theWiseAman]

This is how it will work now. I have removed the nthroot() as suggested and used RootsVector <: AbstractArray. I even have generalised RootsVector for both z::Complex as well as z::Real.

julia> roots = RootsVector(1+0im, 4)
4-element RootsVector{Complex{Int64}}:
                     1.0 + 0.0im
   6.123233995736766e-17 + 1.0im
                    -1.0 + 1.2246467991473532e-16im
 -1.8369701987210297e-16 - 1.0im

julia> roots = RootsVector(1+0im, 4, 'm')
4-element RootsVector{Complex{Int64}}:
                     1.0 + 0.0im
   6.123233995736766e-17 + 1.0im
                    -1.0 + 1.2246467991473532e-16im
 -1.8369701987210297e-16 - 1.0im

julia> roots = RootsVector(1+0im, 4, 's')
1-element RootsVector{Complex{Int64}}:
 1.0 + 0.0im

julia> roots = RootsVector(1, 4)
1-element RootsVector{Int64}:
 1.0

julia> roots = RootsVector(1, 4, 'm')
2-element RootsVector{Int64}:
  1.0
 -1.0

julia> roots = RootsVector(1, 4, 's')
1-element RootsVector{Int64}:
 1.0

julia> roots = RootsVector(1, 4, 'm')
2-element RootsVector{Int64}:
  1.0
 -1.0

julia> roots = RootsVector(1, 4, 's')
1-element RootsVector{Int64}:
 1.0

Need feedback on z::Real && z < 0. I am currently returning NaN currently when n::even

julia> roots = RootsVector(-1, 4, 'm')
1-element RootsVector{Int64}:
 NaN

julia> roots = RootsVector(-1, 4)
1-element RootsVector{Int64}:
 NaN

julia> roots = RootsVector(-1, 4, 's')
1-element RootsVector{Int64}:
 NaN

julia> roots = RootsVector(8, 3, 'm')
1-element RootsVector{Int64}:
 2.0

julia> roots = RootsVector(8, 3, 's')
1-element RootsVector{Int64}:
 2.0

julia> roots = RootsVector(8, 3)
1-element RootsVector{Int64}:
 2.0

julia> roots = RootsVector(-8, 3)
1-element RootsVector{Int64}:
 -2.0

julia> roots = RootsVector(-8, 3, 'm')
1-element RootsVector{Int64}:
 -2.0

julia> roots = RootsVector(-8, 3, 's')
1-element RootsVector{Int64}:
 -2.0

[LilithHafner]

There are a lot of features here. Perhaps it makes sense to prototype in a package first and see if there are many use cases where these features are needed.

[theWiseAman]

sqrt has a branch cut along the negative real line. See also https://en.wikipedia.org/wiki/Branch_point#Branch_cuts

I think for sqrt branch cut is along the imaginary axis.

So, I have implemented single-valued function mode for RootsVector() represented as option = 's'. I am using angle() which returns θ ∈ −π < θ < π, so I think branch cut would be along θ = π/n && θ = 0 or θ = 0 && θ = -π/n. Basically for single mode it is either returning 0 < RootsVector(z, n, 's') < π/n or -π/n < RootsVector(z, n, 's') < 0.

Also RootsVector(z, n, 's') is consistent with sqrt()

julia> x = -1
-1

julia> Δ = 1e-9im
0.0 + 1.0e-9im

julia> RootsVector(x+Δ, 3, 's')
1-element RootsVector{ComplexF64}:
 0.5000000002886753 + 0.8660254036177719im

julia> RootsVector(x-Δ, 3, 's')
1-element RootsVector{ComplexF64}:
 0.5000000002886753 - 0.8660254036177719im

julia> RootsVector(x+Δ, 2, 's')
1-element RootsVector{ComplexF64}:
 5.000001026025254e-10 + 1.0im

julia> RootsVector(x-Δ, 2, 's')
1-element RootsVector{ComplexF64}:
 5.000001026025254e-10 - 1.0im

julia> RootsVector(x+Δ, 2, 's') - RootsVector(x-Δ,2,'s')
1-element Vector{ComplexF64}:
 0.0 + 2.0im

[theWiseAman]

There are a lot of features here. Perhaps it makes sense to prototype in a package first and see if there are many use cases where these features are needed.

@LilithHafner I think the RootsVector() is fully complete and ready to merge. It just needs a good documentation and tests. From my side I have implemented all the tools that might be needed in different real life physical problems.

Suggestions are appreciated. Please confirm if it's good to merge, then I would like to pull the JuliaLang:master branch and make my PR available for merging. Any suggestions on rebasing this branch?

[LilithHafner]

There is an infinite set of possible features to add to the language, but it's nice to keep the language itself small. Just because a feature is implemented tested and documented doesn't mean it is appropriate for Base. For example, which use cases motivate the inclusion of both an 's' and 'm' mode?

Feel free to merge master into this branch or rebase this branch on top of the baster branch. I typically prefer it when PR authors merge so that I still have the option to view "changes since last review", but other folks prefer rebasing. Either way the entire PR will get squashed into a single commit to the Julia master branch if it is merged, so there is no need to be particularly fussy about git history, what matters most is the content of the PR, not its git history.

As an aside, typically folks use backticks (`) for literal code snippets like sqrt and nthroots(4, 2), not for more general concepts like documentation, tests, and rebasing.

[thchr]

Do we have a function which returns roots of unity?

cispi.(range(0, 2-2/n, n))

(if performance matters, it might make sense to do it by repeated multiplication by cispi(2/n): but performance is likely irrelevant here).

I think n roots might be useful, so we should have that option to extract it.

What are some examples of this potential usefulness beyond pedagogy?

[theWiseAman]

@LilithHafner
Resolves issues #47565 and #36534.
Can you merge this if changes are not needed? Else could you suggest any changes?

[theWiseAman]

What are some examples of this potential usefulness beyond pedagogy?

I am Physics undergrad students, so I thought it might be useful. But I don't have the vantage point to comment on it.

Just mention me in the comments if some changes are needed. I would like my first PR to complete all the steps. I would be happy to incorporate the changes the community wants.

[LilithHafner]

There are still several unresolved comments including

1. Where to put the branch cut?
2. Is it good to return multiple results?
3. What are some examples of this potential usefulness beyond pedagogy?

These are more design decisions that need to be made than things that need changing in this PR.

[stevengj]

I think this is a bad idea. The whole point of an nthroot function is to return the real root for real arguments, which exists even for negative reals if you have odd powers, where it differs from the principal root returned by x^(1/n).

I don't think we want a function to return an array (or iterator) of all roots for complex arguments. I don't think nthroot should support complex arguments at all to start with. We can potentially do something along the lines of #36534 (comment) later on, but to begin with it should support only ::Real arguments, which should be less controversial.

If you want to work on this, I would suggest closing this PR and starting over from scratch. In general, if you want to resolve an issue you should implement the functionality discussed in the issue, not some completely different function.