Support DynamicPolynomials.jl

[dlfivefifty]

This supports combining this package with symbolic computation. It's surprisingly fast! 3.5x faster than mathematica in fact:

julia> @time chebyshevt(10_000,x)
  0.078005 seconds (3.84 M allocations: 141.068 MiB, 18.68% gc time)
1 - 50000000x² + 416666650000000x⁴ - 5382805582903628800x⁶ - 2295839942961001984x⁸ - 5881664067181658112x¹⁰ + 9071434676160905216x¹² + 2726131615473205248x¹⁴ - 8311353073181097984x¹⁶ - 7478286491391623168x¹⁸ - 6167377969439309824x²⁰ - 749395590768492544x²² - 3810475745168326656x²⁴ + 5035205751279190016x²⁶ + 5745040218969341952x²⁸ - 5921272124512665600x³⁰ + 3935048435299778560x³² + 5526248795295186944x³⁴ - 4574133298292326400x³⁶ + 379621782652452864x³⁸ - 7095975366782615552x⁴⁰ - 9126263169889599488x⁴² - 5118059501529858048x⁴⁴ + 6924284427082137600x⁴⁶ + 7205196453839372288x⁴⁸ + 1729382256910270464x⁵⁰ + 2882303761517117440x⁵² - 4611686018427387904x⁵⁴ + 6917529027641081856x⁵⁶

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 92.15%. Comparing base (50a37b4) to head (001c4a7).
Report is 1 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main     #214   +/-   ##
=======================================
  Coverage   92.14%   92.15%           
=======================================
  Files          18       19    +1     
  Lines        1858     1860    +2     
=======================================
+ Hits         1712     1714    +2     
  Misses        146      146           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[blegat]

JuMP expressions are also AbstractMutable, I'm wondering if this works now:

using JuMP
model = Model()
@variable(model, x)
chebyshevt(5, x)

I'm getting

julia> chebyshevt(5, x)
ERROR: MethodError: Cannot `convert` an object of type QuadExpr to an object of type AffExpr
The function `convert` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  GenericAffExpr{V, K}(::Any, Pair...) where {K, V, N}
   @ JuMP ~/.julia/packages/JuMP/i68GU/src/aff_expr.jl:205
  GenericAffExpr{V, K}(::Any, ::AbstractArray{<:Pair}) where {K, V}
   @ JuMP ~/.julia/packages/JuMP/i68GU/src/aff_expr.jl:198
  (::Type{GenericAffExpr{CoefType, VarType}} where {CoefType, VarType})(::Any, ::Any)
   @ JuMP ~/.julia/packages/JuMP/i68GU/src/aff_expr.jl:140
  ...

Stacktrace:
 [1] initiateforwardrecurrence(N::Int64, A::LazyArrays.ApplyArray{…}, B::FillArrays.Zeros{…}, C::FillArrays.Ones{…}, x::VariableRef, μ::AffExpr)
   @ RecurrenceRelationshipArrays ~/.julia/packages/RecurrenceRelationshipArrays/SCsob/src/recurrence.jl:30
 [2] unsafe_getindex
   @ ~/.julia/dev/ClassicalOrthogonalPolynomials/ext/ClassicalOrthogonalPolynomialsMutableArithmeticsExt.jl:5 [inlined]
 [3] chebyshevt(n::Int64, z::VariableRef)
   @ ClassicalOrthogonalPolynomials ~/.julia/dev/ClassicalOrthogonalPolynomials/src/classical/chebyshev.jl:64
 [4] top-level scope
   @ REPL[9]:1
Some type information was truncated. Use `show(err)` to see complete types.

The issue is that the types are tricky to get right, you first have a JuMP.VariableRef, then a JuMP.AffExpr then a JuMP.QuadExpr then it stabilises at JuMP.NonlinearExpr. So you would need to run a fixed point with MutableArithmetics.promote_operation ^^

[dlfivefifty]

The offending line is here:

https://github.com/JuliaApproximation/RecurrenceRelationshipArrays.jl/blob/bac197ddaa3bd5a1fbdb0f838d98e8658d21a734/src/recurrence.jl#L30

You can see it makes the mistake of assuming the output type is just the promotion of the input types. Probably a better design might be:

    p1 = muladd(A[1],x,B[1])*μ
    T = typeof(p1)
     p0 = convert(T, μ)

Though some thought would be needed for the N == 0 special case...

[blegat]

You probably don't want the type of the result to depend on N so it's fine if you use the type corresponding to a larger N