Poles and roots calculation failure

[Fabi293]

using Polynomials

#setprecision(BigFloat,12)  # also wrong for small precisions

case = 1
if case == 1  # UndefRefError
	numcoeffs = Complex{BigFloat}[-0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, -6.4095e+257 - 3.314e+243im, -5.6118e+244 - 3.9514e+230im, -4.0102e+236 - 2.0735e+222im, -3.0725e+223 - 2.1631e+209im, -1.0975e+215 - 5.6751e+200im, -7.2072e+201 - 5.0751e+187im, -1.7167e+193 - 8.877e+178im, -9.3936e+179 - 6.6134e+165im, -1.678e+171 - 8.6782e+156im, -7.3448e+157 - 5.1693e+143im, -1.0499e+149 - 5.43e+134im, -3.4468e+135 - 2.4259e+121im, -4.1052e+126 - 2.1231e+112im, -8.9846e+112 - 6.3225e+98im, -9.1725e+103 - 4.742e+89im, -1.0036e+90 - 7.0656e+75im, -8.9646e+80 - 4.6354e+66im]
	dencoeffs = Complex{BigFloat}[0.0 + 0.0im, 0.0 + 0.0im, -0.0 + 0.0im, -0.0 + 0.0im, 0.0 + 0.0im, 0.0 + 0.0im, -3.0837e+258 + 0.0im, -3.5997e+245 - 9.309e+230im, -1.9292e+237 - 7.0217e+217im, -1.9709e+224 - 5.0951e+209im, -5.2803e+215 - 3.2932e+196im, -4.6214e+202 - 1.1952e+188im, -8.2587e+193 - 6.4383e+174im, -6.0237e+180 - 1.5577e+166im, -8.0763e+171 - 6.7125e+152im, -4.711e+158 - 1.2183e+144im, -5.0507e+149 - 3.937e+130im, -2.2106e+136 - 5.7173e+121im, -1.9752e+127 - 1.2316e+108im, -5.7628e+113 - 1.4904e+99im, -4.4151e+104 - 1.6057e+85im, -6.4393e+90 - 1.6651e+76im, -4.3167e+81 + 0.0im]
else  # MethodError
	numcoeffs_ = collect(1:20)
	numcoeffs = Complex{BigFloat}[]
	numcoeffs = Complex{BigFloat}.(numcoeffs_)

	dencoeffs_ = collect(21:40)
	dencoeffs = Complex{BigFloat}[]
	dencoeffs = Complex{BigFloat}.(dencoeffs_)
end

p = Polynomial(numcoeffs)
q = Polynomial(dencoeffs)
r = p//q

poles(r)
#roots(r)

If the given code with case = 1 is executed, I get the error

ERROR: UndefRefError: access to undefined reference
Stacktrace:
[1] getindex
@ .\essentials.jl:917 [inlined]
[2] getindex
@ .\array.jl:930 [inlined]
[3] getindex
@ USER.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\triangular.jl:286 [inlined]
[4] copyto!(A::LinearAlgebra.UpperTriangular{…}, B::LinearAlgebra.UpperTriangular{…})
@ LinearAlgebra USER.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\triangular.jl:551
[5] copy
@ USER.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\triangular.jl:190 [inlined]
[6] refine_uvw!(u::Polynomials.MutableDenseViewPolynomial{…}, v::Polynomials.MutableDenseViewPolynomial{…}, w::Polynomials.MutableDenseViewPolynomial{…}, p::Polynomials.MutableDenseViewPolynomial{…}, q::Polynomials.MutableDenseViewPolynomial{…}, uv::Polynomials.MutableDenseViewPolynomial{…}, uw::Polynomials.MutableDenseViewPolynomial{…})
@ Polynomials.NGCD USER.julia\packages\Polynomials\6i39P\src\polynomials\ngcd.jl:543
[7] ngcd(p::Polynomials.MutableDenseViewPolynomial{…}, q::Polynomials.MutableDenseViewPolynomial{…}; atol::BigFloat, rtol::BigFloat, satol::BigFloat, srtol::BigFloat, scale::Bool, λ::BigFloat, minⱼ::Int64)
@ Polynomials.NGCD USER.julia\packages\Polynomials\6i39P\src\polynomials\ngcd.jl:295
[8] ngcd(p::Polynomials.MutableDenseViewPolynomial{…}, q::Polynomials.MutableDenseViewPolynomial{…})
@ Polynomials.NGCD USER.julia\packages\Polynomials\6i39P\src\polynomials\ngcd.jl:228
[9] ngcd(::Polynomial{Complex{BigFloat}, :x}, ::Polynomial{Complex{BigFloat}, :x}; kwargs::@kwargs{})
@ Polynomials USER.julia\packages\Polynomials\6i39P\src\polynomials\ngcd.jl:51
[10] ngcd
@ USER.julia\packages\Polynomials\6i39P\src\polynomials\ngcd.jl:11 [inlined]
[11] #uvw#130
@ USER.julia\packages\Polynomials\6i39P\src\polynomials\standard-basis\standard-basis.jl:349 [inlined]
[12] uvw
@ USER.julia\packages\Polynomials\6i39P\src\polynomials\standard-basis\standard-basis.jl:348 [inlined]
[13] uvw(p::Polynomial{Complex{BigFloat}, :x}, q::Polynomial{Complex{BigFloat}, :x}; method::Symbol, kwargs::@kwargs{})
@ Polynomials USER.julia\packages\Polynomials\6i39P\src\polynomials\standard-basis\standard-basis.jl:343
[14] uvw
@ USER.julia\packages\Polynomials\6i39P\src\polynomials\standard-basis\standard-basis.jl:343 [inlined]
[15] lowest_terms(pq::RationalFunction{Complex{BigFloat}, :x, Polynomial{Complex{…}, :x}}; method::Symbol, kwargs::@kwargs{})
@ Polynomials USER.julia\packages\Polynomials\6i39P\src\rational-functions\common.jl:424
[16] lowest_terms
@ USER.julia\packages\Polynomials\6i39P\src\rational-functions\common.jl:420 [inlined]
[17] #poles#184
@ USER.julia\packages\Polynomials\6i39P\src\rational-functions\common.jl:436 [inlined]
[18] poles(pq::RationalFunction{Complex{BigFloat}, :x, Polynomial{Complex{BigFloat}, :x}})
@ Polynomials USER.julia\packages\Polynomials\6i39P\src\rational-functions\common.jl:435
[19] top-level scope
@ PATH\juliaBugs\pol2.jl:23
Some type information was truncated. Use show(err) to see complete types.

While trying to make case 1 simpler, I got a different error (case = 2):

ERROR: MethodError: no method matching eigvals!(::Matrix{Complex{BigFloat}})
The function eigvals! exists, but no method is defined for this combination of argument types.

Closest candidates are:
eigvals!(::AbstractMatrix{T}, ::LinearAlgebra.Cholesky{T}; sortby) where T<:Number
@ LinearAlgebra USER.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\symmetriceigen.jl:260
eigvals!(::StridedMatrix{T}, ::LinearAlgebra.LU{T, <:StridedMatrix{T} where T}; sortby) where T
@ LinearAlgebra USER.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\symmetriceigen.jl:286
eigvals!(::LinearAlgebra.Hermitian{T, S}, ::LinearAlgebra.Hermitian{T, S}; sortby) where {T<:Union{ComplexF64, ComplexF32}, S<:(StridedMatrix{T} where T)}
@ LinearAlgebra USER.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\symmetriceigen.jl:246
...

Stacktrace:
[1] eigvals(A::Matrix{Complex{BigFloat}}; kws::@kwargs{})
@ LinearAlgebra USER.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\eigen.jl:343
[2] eigvals(A::Matrix{Complex{BigFloat}})
@ LinearAlgebra USER.julia\juliaup\julia-1.11.3+0.x64.w64.mingw32\share\julia\stdlib\v1.11\LinearAlgebra\src\eigen.jl:343
[3] roots(p::Polynomials.MutableDenseViewPolynomial{Polynomials.StandardBasis, Complex{BigFloat}, :x}; kwargs::@kwargs{})
@ Polynomials USER.julia\packages\Polynomials\6i39P\src\polynomials\standard-basis\standard-basis.jl:499
[4] roots
@ USER.julia\packages\Polynomials\6i39P\src\polynomials\standard-basis\standard-basis.jl:480 [inlined]
[5] pejorative_manifold(p::Polynomial{Complex{…}, :x}; method::Symbol, θ::Float64, ρ::Float64, ϕ::Int64, kwargs::@kwargs{})
@ Polynomials.Multroot USER.julia\packages\Polynomials\6i39P\src\polynomials\multroot.jl:160
[6] pejorative_manifold
@ USER.julia\packages\Polynomials\6i39P\src\polynomials\multroot.jl:140 [inlined]
[7] multroot(p::Polynomial{Complex{BigFloat}, :x}; verbose::Bool, kwargs::@kwargs{})
@ Polynomials.Multroot USER.julia\packages\Polynomials\6i39P\src\polynomials\multroot.jl:114
[8] multroot
@ USER.julia\packages\Polynomials\6i39P\src\polynomials\multroot.jl:101 [inlined]
[9] #poles#184
@ USER.julia\packages\Polynomials\6i39P\src\rational-functions\common.jl:438 [inlined]
[10] poles(pq::RationalFunction{Complex{BigFloat}, :x, Polynomial{Complex{BigFloat}, :x}})
@ Polynomials USER.julia\packages\Polynomials\6i39P\src\rational-functions\common.jl:435
[11] top-level scope
@ PATH\juliaBugs\pol2.jl:23
Some type information was truncated. Use show(err) to see complete types.

The same error comes when calculating the roots.

[jverzani]

Hmm, for these the issue seems to lie in the call to multroot in the poles function. This defaults to a ":direct" method, which is failing with these big coefficients, whereas the alternative ":iterative" method works. However, there is currently no means to specify that through poles. The simple patch is to allow a keyword argument to slide through, but that isn't so helpful unless an error message suggests this implementation detail. I'll add this first and will have to think about the latter.

[Fabi293]

Thanks for the fix! But similar errors still exist for the roots method.

[jverzani]

Ohh, thanks for pointing this out. I'll try and get to it soon.