fix evalpoly type instability and 0-length case

[stevengj]

A few fixes/improvements to the evalpoly function:

1. Fixes evalpoly(x,()) gives unhelpful error. #56699: returns zero for an empty tuple or array of coefficients
2. Fixes a type-instability in evalpoly(x, ::Vector)
3. Adds a missing @inbounds annotation
4. Trivial: Renames a local variable ex that was-copy-pasted from the macro version (where ex referred to an expression) but makes no sense here for the numeric sum.
5. Removes some extraneous methods (by tightening a method signature, as suggested by @nsajko) — should still dispatch to the same implementations as before.

Probably should be backported.

[nsajko]

I'm not sure if opening this PR was polite, considering a newcomer submitted #56699 trying to fix the same issue, and the issue is marked as "good first issue". In any case, that PR is now closed, so here are some comments.

[stevengj]

(I hadn't noticed the other PR until after I submitted this one, unfortunately.)

[stevengj]

CI failures look unrelated.

[stevengj]

Yay, green. Good to merge?

[oscardssmith]

why is triage added here?

[stevengj]

Because it's been ready to merge for 2 months and needs someone to review it? Is there a better label?

[LilithHafner]

This makes the existing use of @generated on line 95 invalid because the two branches now return different types:

julia> evalpoly(1.0, (1,))
1

julia> Base.Math._evalpoly(1.0, (1,))
1.0

julia> @less evalpoly(1, (1,))
function evalpoly(x, p::Tuple)
    if @generated
        N = length(p.parameters::Core.SimpleVector)
        ex = :(p[end])
        for i in N-1:-1:1
            ex = :(muladd(x, $ex, p[$i]))
        end
        ex
    else
        _evalpoly(x, p)
    end
end

[stevengj]

Oh ye gods, børked rebase.

[giordano]

Git history is messed up now.

[stevengj]

Should be fixed now.

[stevengj]

@LilithHafner, I fixed the issue with the @generated branches you noted in #56707 (review), and added a test.

[stevengj]

Currently getting:

 Test threw exception
  Expression: #= /cache/build/tester-amdci4-13/julialang/julia-master/julia-f7b9e88614/share/julia/stdlib/v1.12/LinearAlgebra/test/matmul.jl:647 =# @evalpoly(A33, 1.0I, 1.0I) == I + A33
  MethodError: Cannot `convert` an object of type LinearAlgebra.UniformScaling{Float64} to an object of type Matrix{ComplexF64}
  The function `convert` exists, but no method is defined for this combination of argument types.

This could be fixed by replacing oftype(one(x) * p[end], p[end]) with simply one(x) * p[end]. However, in that case I'm concerned that JuliaLang/LinearAlgebra.jl#1161 will get worse — it will then start silently giving the wrong answer for matrix x and scalar coefficients, because of the weird behavior of muladd (JuliaLang/LinearAlgebra.jl#1162) in this case.

So maybe this PR is blocked by JuliaLang/LinearAlgebra.jl#1161 being fixed within LinearAlgebra.jl?

[LilithHafner]

This now breaks the identity evalpoly(x, (y,)) === y. e.g:

julia> evalpoly(:sploo, (:bork,))
:bork # now throws

julia> evalpoly(1.0, (true,))::Bool
true # now 1.0

julia> evalpoly(rand(2,2), (1,))
1 # This should possibly be Float64[1 0; 0 1] now throws

According to the docstring, evalpoly(x, (y,)) should return x^0 * y which, after this PR, in the first two cases it does and in the third case it hits JuliaLang/LinearAlgebra.jl#1161.

I think that "breakage" is a bugfix.

[stevengj]

This now breaks the identity evalpoly(x, (y,)) === y

Not sure why this identity should be expected from the definition?

I think that "breakage" is a bugfix.

I guess you agree with the new behavior?

[LilithHafner]

Not sure why this identity should be expected from the definition?

If you write out a polynomial in the form ax^2+bx+c, rather than a_2x^2+a_1x^1+a_0x^0, then the degree 0 case is just a which would imply the result would be egal to a.

I guess you agree with the new behavior?

Yes, I do. Regardless of the above.

[stevengj]

If you write out a polynomial in the form ax^2+bx+c

This is not the definition that is given in the docs for evalpoly. It is defined as $\sum_k x^{k-1} p[k]$.

(Note also that the definition right-multiplies p[k]. I'm not sure if there is someone out there who prefers left-multiplied p[k] for non-commutative algebras, e.g. matrix x and p[k]? They'll have to define their own evalpoly, I guess.)

[fingolfin]

There were CI failures, e.g.

mith-C07ZM05NJYVY-0/julialang/julia-master/julia-f7b9e88614/share/julia/stdlib/v1.12/LinearAlgebra/test/matmul.jl:647
  Test threw exception
  Expression: #= /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-grannysmith-C07ZM05NJYVY.0/build/default-grannysmith
-C07ZM05NJYVY-0/julialang/julia-master/julia-f7b9e88614/share/julia/stdlib/v1.12/LinearAlgebra/test/matmul.jl:647 =# @evalpoly(A33, 1.0I, 1.0I) == I + A33
  MethodError: Cannot `convert` an object of type LinearAlgebra.UniformScaling{Float64} to an object of type Matrix{ComplexF64}
  The function `convert` exists, but no method is defined for this combination of argument types.
  Closest candidates are:
    Matrix{T}(::LinearAlgebra.UniformScaling, !Matched::Tuple{Int64, Int64}) where T
     @ LinearAlgebra /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-grannysmith-C07ZM05NJYVY.0/build/default-grannys
mith-C07ZM05NJYVY-0/julialang/julia-master/julia-f7b9e88614/share/julia/stdlib/v1.12/LinearAlgebra/src/uniformscaling.jl:423
    Matrix{T}(::LinearAlgebra.UniformScaling, !Matched::Integer, !Matched::Integer) where T
     @ LinearAlgebra /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-grannysmith-C07ZM05NJYVY.0/build/default-grannys
mith-C07ZM05NJYVY-0/julialang/julia-master/julia-f7b9e88614/share/julia/stdlib/v1.12/LinearAlgebra/src/uniformscaling.jl:431
    convert(::Type{T}, !Matched::AbstractArray) where T<:Array
     @ Base array.jl:614

[stevengj]

Blocked by JuliaLang/LinearAlgebra.jl#1161 as noted above.

[MasonProtter]

I worry a bit about this causing performance regressions (especially in cases where multiplication is not cheap), but I also get that it's likely important.

Ideally we'd be able to do something in the type domain here rather than an actual multiplication, but oh well i guess. This is probably less brittle, and won't regress in most cases.

[stevengj]

Hopefully the compiler can optimize out one(x) * p[end] in most cases, and just compute the type, since the result is not used.

[fingolfin]

Only if the multiplication is simple enough and free of side effects like memory allocation? If the scalar is eg a BigInt it may not be able to optimize this away?