improved implementations of composition of polynomials

[nsajko]

1. it seems that algorithms with less computational complexity than mere Horner's scheme exist, for example (covers some background in the introduction): https://ens-lyon.hal.science/ensl-00546102/document
2. There should be special code for the case where one or both polynomials have their (small) degree in the type domain (ImmutablePolynomial)

[jverzani]

Thanks! I'll have a look when I get a chance.

[nsajko]

Regarding my point two above, here's a proof-of-concept implementation for composing a Polynomial with a degree one ImmutablePolynomial:

function weighted_horner(x, p, f)
  deg = length(p) - 1
  w = f(Val(:init), deg)
  out = (p[end] * w) * one(x)
  for i ∈ (deg - 1):-1:0
    w *= f(Val(:mul), i)
    w = div(w, f(Val(:div), i), RoundToZero)
    out = muladd(out, x, p[begin + i] * w)
  end
  out
end

struct PolynomialCompositionHelper
  k::Int
end
bin_coef_bottom(h::PolynomialCompositionHelper) = h.k
function (h::PolynomialCompositionHelper)(::Val{:init}, deg::Int)
  k = bin_coef_bottom(h)
  binomial(deg + k, k)
end
(::PolynomialCompositionHelper)(::Val{:mul}, deg::Int) = deg + 1
(h::PolynomialCompositionHelper)(::Val{:div}, deg::Int) =
  deg + 1 + bin_coef_bottom(h)

function polynomial_composed_with_shifted_identity(
  p::Vector{T},
  shift::T,
) where {T}
  len = length(p)
  out = Vector{T}(undef, len)
  for k ∈ 0:(len - 1)
    h = PolynomialCompositionHelper(k)
    out[begin + k] = weighted_horner(shift, (@view p[(begin + k):end]), h)
  end
  out
end

function compose_polynomial_with_linear!(p::Vector{T}, scale::T) where {T}
  s = scale
  for i ∈ 1:(length(p) - 1)
    p[begin + i] *= s
    s *= scale
  end
  p
end

# Composition `p ∘ r`. `p` and `r` are the coefficients of polynomials
# in the standard basis.
composed_polynomials(p::Vector{T}, r::NTuple{2,T}) where {T} =
  # Why this works:
  #
  # 1. `r` can be decomposed like so: `r = s ∘ l`, where
  #    `r(x) = a + b*x`, `s(x) = a + x`, `l(x) = b*x`
  #
  # 2. `p ∘ r = p ∘ (s ∘ l) = (p ∘ s) ∘ l`
  compose_polynomial_with_linear!(
    polynomial_composed_with_shifted_identity(p, first(r)),
    last(r),
  )

using Polynomials

(p::Polynomial{T, v})(r::ImmutablePolynomial{T, v, 2}) where {T, v} =
  Polynomial{T, v}(composed_polynomials(coeffs(p), coeffs(r)))

I didn't compare accuracy yet, however the results appear to be approximately correct:

julia> include("polynomial_composition.jl")  # includes the above code

julia> p = Polynomial(rand(6));

julia> q = ImmutablePolynomial((rand(), rand()))
ImmutablePolynomial(0.7923048129450314 + 0.4495999592340951*x)

julia> p(q)  # new implementation
Polynomial(2.2449893995110872 + 2.4472444477738295*x + 1.902292181529416*x^2 + 0.8098154846551124*x^3 + 0.18248093200464596*x^4 + 0.017903450489708504*x^5)

julia> p(Polynomial(q))  # old implementation
Polynomial(2.244989399511087 + 2.447244447773829*x + 1.9022921815294158*x^2 + 0.8098154846551123*x^3 + 0.18248093200464593*x^4 + 0.017903450489708504*x^5)

The above implementation is faster than what's currently available in Polynomials.jl, however I'll wait until your upcoming PR is merged before making concrete comparisons.

[jverzani]

Thanks! I'll have to look this over to see how to phase it in. It seems like there should be some big performance optimizations when the immutable polynomial type is involved.

[nsajko]

After playing some more with my implementation above, I realized that the approach is possibly problematic: in weighted_horner, the weight w is mathematically an integer, however it's huge for large-degree polynomials, and f(Val(:init), deg) overflows Int starting when deg is somewhere between 60 and 70.

So I guess this should be adapted slightly and restricted to something like Polynomial{<:Union{AbstractFloat,BigInt}}.