Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix evalpoly type instability and 0-length case #56707

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open

fix evalpoly type instability and 0-length case #56707

wants to merge 10 commits into from

Conversation

stevengj
Copy link
Member

@stevengj stevengj commented Nov 28, 2024 •
edited
Loading

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.

@stevengj stevengj added maths Mathematical functions bugfix This change fixes an existing bug labels Nov 28, 2024
@Sagar-Ramchandani Sagar-Ramchandani mentioned this pull request Nov 28, 2024
Closed
LilithHafner
LilithHafner reviewed Nov 28, 2024
View reviewed changes
base/math.jl Outdated Show resolved Hide resolved
nsajko
nsajko suggested changes Nov 29, 2024
View reviewed changes
Copy link
Contributor

@nsajko nsajko left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

base/math.jl Outdated Show resolved Hide resolved
base/math.jl Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
base/math.jl Outdated Show resolved Hide resolved
base/math.jl Show resolved Hide resolved
base/math.jl Show resolved Hide resolved
@stevengj
Copy link
Member Author

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

@stevengj
Copy link
Member Author

stevengj commented Dec 2, 2024

CI failures look unrelated.

@stevengj
Copy link
Member Author

Yay, green. Good to merge?

@stevengj stevengj added the triage This should be discussed on a triage call label Dec 31, 2024
@oscardssmith
Copy link
Member

why is triage added here?

@stevengj
Copy link
Member Author

stevengj commented Jan 1, 2025

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

@stevengj stevengj removed the triage This should be discussed on a triage call label Jan 1, 2025
LilithHafner
LilithHafner reviewed Jan 1, 2025
View reviewed changes
Copy link
Member

@LilithHafner LilithHafner left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Member Author

stevengj commented Jan 2, 2025

Oh ye gods, børked rebase.

@giordano
Copy link
Contributor

giordano commented Jan 2, 2025

Git history is messed up now.

@stevengj stevengj force-pushed the sgj/evalpoly_fixes branch from c6852d2 to 05cfd26 Compare January 2, 2025 15:13
@stevengj
Copy link
Member Author

stevengj commented Jan 2, 2025

Should be fixed now.

@stevengj
Copy link
Member Author

stevengj commented Jan 2, 2025

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

@stevengj
Copy link
Member Author

stevengj commented Jan 2, 2025 •
edited
Loading

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
Copy link
Member

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
Copy link
Member Author

stevengj commented Jan 3, 2025 •
edited
Loading

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
Copy link
Member

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
Copy link
Member Author

stevengj commented Jan 3, 2025 •
edited
Loading

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.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@LilithHafner LilithHafner LilithHafner left review comments

@nsajko nsajko nsajko requested changes

Assignees
No one assigned
Labels
bugfix This change fixes an existing bug maths Mathematical functions
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.

evalpoly(x,()) gives unhelpful error.
5 participants
@stevengj @oscardssmith @giordano @LilithHafner @nsajko