## Connection and Linearization

## A connection between polynomial systems `P` and `Q` is a function `C(m,n)` satisfying
## p_n = ∑C(m,n) q_m or  pᵢ = ∑ Cʲᵢqⱼ (if unicode is helpful)
##
## From this a polynomial in P = ∑ aᵢ⋅pᵢ can be represented in a Q through
##
##  ∑ᵢ aᵢ⋅pᵢ = ∑ᵢ aᵢ  ∑ⱼ Cʲᵢqⱼ  = ∑ᵢ (∑ⱼ (aᵢ⋅Cʲᵢ)) qⱼ or ∑ⱼ (∑ᵢ aᵢCʲᵢ) qⱼ
##
## The choice of which depends on how the connection coefficients are presented.
##
## A connection between P and Q defines a `convert` method `Base.convert(Q,  p::P)`
## Conversion can also be achieved through polynomial evaluation. For example
## if `x=variable(Q)` then `Base.convert(Q,p) = p(x)` will also convert `p` into
## the basis Q. This is helpful when Q is a standard basis polynomial.
##
## However, if Q is not a standard basis polynomial, it isn't helpful unless
## one can multiply monomomials within Q. For standard basis polynomials this
## is known, as  xⁱ⋅xʲ = xⁱ⁺ʲ.
##
## For general Q though, one needs a linearization: a function α with
##
##  qᵢqⱼ = ∑_k α(k,i,j) q_k
##
## with a given α, multiplication can be defined. The linearization can be defined
## through composition: If P is the standard basis, connect Q to P, multiply ini P, then connect
## to Q.
##
## In short, if there is a linearization then multiplication can be defined and conversion defined through
## polynomial evaluation
## If there is a connection defined then conversion is defined and multiplication can be defined through composition
##


## For `AbstractCCOP` polynomials, Koepf and Schmersau present recursive formulas for the connection coefficients
## between Q and the standard basis, the standard and Q, and between P and Q when σ and σ̂ are the same. (This is the case when a polynomial system is parameterized, like Laguerre, Gegenbauer, Bessel, and  Jacobi.)

ConvertibleTypes = Union{AbstractCOP, Polynomials.StandardBasisPolynomial,  FallingFactorial}
kn(::Type{P},  n) where{P <: Union{Polynomials.StandardBasisPolynomial,  FallingFactorial}} = one(eltype(one(P)))
k1k0(::Type{P},  n) where{P <: Union{Polynomials.StandardBasisPolynomial,  FallingFactorial}} = one(eltype(one(P)))
knk_1(::Type{P},  n) where{P <: Union{Polynomials.StandardBasisPolynomial,  FallingFactorial}} = one(eltype(one(P)))



## If  the connection  coefficients for  (P,Q) satisfy c₀⋅ C̃ⁱⱼ + c₁⋅C̃ⁱ⁺¹ⱼ  + c₂⋅*C̃ⁱ⁺²ⱼ = 0 for some computable
##  values c₀,c₁,c₂ then this  will implement  `convert(Q,p::P)`
function _convert_cop(::Type{Q}, p::P) where {P <: ConvertibleTypes,
                                               Q <: ConvertibleTypes}

    d = degree(p)
    T,S = eltype(Q), eltype(p)  #
    R = typeof(one(T) * p[0]/1)
    as = zeros(R, 1+d)

    λⱼⱼ = one(R)  * k0(P) / k0(Q)


    for j in  0:d

        λ = λⱼⱼ

        λⱼⱼ *= k1k0(P,j) / k1k0(Q,j)

        pⱼ = p[j]
        iszero(pⱼ) && continue

        C̃ʲⱼ = one(R)

        as[1+j] += pⱼ * (λ * C̃ʲⱼ)

        C̃ⁱ⁺²ⱼ, C̃ⁱ⁺¹ⱼ = zero(R), C̃ʲⱼ  # for recursion starting  with C̃ʲ⁺¹ⱼ, C̃ʲⱼ
        for i in j-1:-1:0

            λ *= k1k0(Q,i)

            # c₀⋅ C̃ⁱⱼ + c₁⋅C̃ⁱ⁺¹ⱼ  + c₂⋅*C̃ⁱ⁺²ⱼ) = 0
            c₀,c₁,c₂ = connection_m(P, Q, i, j)
            C̃ⁱⱼ = iszero(c₀) ? zero(R) : -one(R)*(c₁*C̃ⁱ⁺¹ⱼ + c₂*C̃ⁱ⁺²ⱼ) / c₀
            as[1+i] += pⱼ * (λ * C̃ⁱⱼ)

            C̃ⁱ⁺²ⱼ, C̃ⁱ⁺¹ⱼ = C̃ⁱ⁺¹ⱼ, C̃ⁱⱼ
        end
    end
    X = Polynomials.indeterminate(Q,p)
    ⟒(Q){R,X}(as)
end

# ## Use FastTransform for conversion, as  possible
# ## FastTransforms.kind2string.(0:15)
# # Legendre<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     T <: AbstractFloat,
#     Q <: Union{Legendre, OrthonormalLegendre},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}}
#     ps = coeffs(p)
#     Q(cheb2leg(ps, normcheb=isorthonormal(P), normleg=isorthonormal(Q)))
# end

# # Chebyshev<--Legendre
# function _convert(::Type{Q}, p::P) where {
#     T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Legendre{T}, OrthonormalLegendre{T}}
# }
#     ps =  coeffs(p)
#     Q(leg2cheb(ps, normleg=isorthonormal(P), normcheb=isorthonormal(Q)) )
# end

# # ultraspherical<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     β, α,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{β}, OrthonormalGegenbauer{β}},
#     P <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}}
# }

#     ps =  coeffs(p)
#     Q( ultra2ultra(ps, α, β, norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end

# # Jacobi<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     γ, δ, α, β,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Jacobi{α, β, T}, OrthonormalJacobi{α, β, T}}
# }

#     ps =  coeffs(p)
#     Q( jac2jac(ps, α, β, γ, δ,  norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end

# # Laguerre<--Laguerre
# function _convert(::Type{Q}, p::P) where {
#     α, β,  T <: AbstractFloat,
#     Q <: Union{Laguerre{β}, OrthonormalLaguerre{β}},
#     P <: Union{Laguerre{α, T}, OrthonormalLaguerre{α, T}}
# }

#     ps =  coeffs(p)
#     Q( lag2lag(ps, α, β,  norm1=isorthonormal(P), norm2=isorthonormal(Q)) )

# end


# # Jacobi<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     γ, δ, α,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}}
# }

#     ps =  coeffs(p)
#     Q( ultra2jac(ps, α, γ, δ,  normultra=isorthonormal(P), normjac=isorthonormal(Q)) )

# end

# # ultraspherical<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     α, γ, δ,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{α, T}, OrthonormalGegenbauer{α, T}},
#     P <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}}
# }

#     ps =  coeffs(p)
#     Q( jac2ultra(ps,γ, δ, α, normjac=isorthonormal(P), normultra=isorthonormal(Q)) )

# end

# # Jacobi<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     γ, δ,  T <: AbstractFloat,
#     Q <: Union{Jacobi{γ, δ}, OrthonormalJacobi{γ, δ}},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}
# }

#     ps =  coeffs(p)
#     Q( cheb2jac(ps,γ, δ, normcheb=isorthonormal(P), normjac=isorthonormal(Q)) )

# end

# # Chebyshev<--Jacobi
# function _convert(::Type{Q}, p::P) where {
#     γ, δ,  T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Jacobi{γ, δ,T}, OrthonormalJacobi{γ, δ,T}}
# }

#     ps =  coeffs(p)
#     Q( jac2cheb(ps,γ, δ, normjac = isorthonormal(P), normcheb=isorthonormal(Q)) )

# end

# # ultraspherical<--Chebyshev
# function _convert(::Type{Q}, p::P) where {
#     α,  T <: AbstractFloat,
#     Q <: Union{Gegenbauer{α}, OrthonormalGegenbauer{α}},
#     P <: Union{Chebyshev{T}, OrthonormalChebyshev{T}}
# }

#     ps =  coeffs(p)
#     Q( cheb2ultra(ps,α, normcheb=isorthonormal(P), normultra=isorthonormal(Q)) )

# end

# # Chebyshev<--ultraspherical
# function _convert(::Type{Q}, p::P) where {
#     α,  T <: AbstractFloat,
#     Q <: Union{Chebyshev, OrthonormalChebyshev},
#     P <: Union{Gegenbauer{α,T}, OrthonormalGegenbauer{α,T}}
# }
#     ps =  coeffs(p)
#     Q( ultra2cheb(ps,α, normultra=isorthonormal(P), normcheb=isorthonormal(Q)) )

# end


##  Koepf and Schmersa Thm 2, p11 connection for P to Q when σ = σ̂
function connection_m(::Type{P},::Type{Q},m,n) where {
    P <: AbstractCCOP,
    Q <: AbstractCCOP}

    a,b,c,d,e = abcde(P)
    ā,b̄,c̄,d̄,ē = abcde(Q)

    (a == ā && b == b̄ && c == c̄) || throw(ArgumentError("This connection assume σ = σ̄"))
    a²,b²,c²,d²,d̄²,e²,ē²,m²,n² = a^2, b^2, c^2, d^2, d̄^2, e^2, ē^2, m^2, n^2

    c0 = -(m-n) ⋅ (a⋅m + d - a + a⋅n) ⋅ (d̄ + 2a⋅m) ⋅ (d̄ + a + 2a⋅m) ⋅ (d̄ + 3a + 2a⋅m)
    c0 *= (d̄ + 2a⋅m  + 2a)^2

    c1 = -d⋅b⋅n⋅d̄ + 2d⋅a⋅m²⋅b + d⋅b⋅d̄ + 2d⋅a⋅m⋅b + 2d⋅ē⋅n⋅a
    c1 += d⋅d̄⋅ē + 2d⋅d̄⋅b⋅m - m⋅b⋅d̄² - e⋅d̄² - 4a²⋅m²⋅e - m²⋅a⋅b⋅d̄ + b⋅n⋅d̄⋅a - 2e⋅d̄⋅a
    c1 += -4a²⋅m⋅e - 4e⋅d̄⋅a⋅m + 2m²⋅a²⋅ē + 2ē⋅a²⋅n² - 2ē⋅a²⋅n - m⋅a⋅b⋅d̄ + 2m⋅d̄⋅ē⋅a
    c1 += 2m⋅ē⋅a² - b⋅n²⋅d̄⋅a
    c1 *= (d̄ + 2a⋅m + 2a) ⋅ (m + 1) ⋅  (d̄ + a + 2a⋅m) ⋅ (d̄ + 3a + 2a⋅m)

    c2 =  -(d̄ + 2a⋅m) ⋅ (m + 1) ⋅ (-a⋅m - 2a + a⋅n - d̄ + d) ⋅ (a⋅m + a⋅n + a + d̄)
    c2a = a⋅b²⋅m² - 4a²⋅m²⋅c - 8a²⋅m⋅c + 2a⋅m⋅b² - 4a⋅d̄⋅m⋅c + m⋅b²⋅d̄ - 4a⋅d̄⋅c - a⋅ē² + a⋅b² - c⋅d̄²
    c2a += b⋅ē⋅d̄ - 4⋅a²⋅c + b²⋅d̄
    c2 *= c2a ⋅ (m + 2)


    return (c0, c1, c2)

end

##  Koepf and Schmersa Thm 4, p14 connection for P, Q=x^n
function connection_m(::Type{P},::Type{Q},m,n) where {
    P <: AbstractCCOP,
    Q <: Polynomials.StandardBasisPolynomial}

    a,b,c,d,e = abcde(P)

    c0 = (m-n) ⋅ (a⋅n + d - a + a⋅m)

    c1 = (m + 1) ⋅ (b⋅m + e)

    c2 = c ⋅ (m + 1) ⋅ (m + 2)

    return (c0, c1, c2)

end

##  Koepf and Schmersa Thm 5, p16 connection for P=x^n, Q
function connection_m(::Type{P},::Type{Q},m,n) where {
    P <: Polynomials.StandardBasisPolynomial,
    Q <: AbstractCCOP}

    ā,b̄,c̄,d̄,ē = abcde(Q)
    ā²,b̄²,c̄²,d̄²,ē²,m²,n² = ā^2, b̄^2, c̄^2, d̄^2, ē^2, m^2, n^2

    c0 =  (n - m) ⋅ (d̄ + 2ā⋅m) ⋅ (d̄ + 3ā + 2ā⋅m) ⋅ (d̄ + ā + 2ā⋅m) ⋅ (d̄ + 2ā⋅m + 2ā)^2

    c1 =  (d̄⋅ē + b̄⋅d̄ +2d̄⋅b̄⋅m + 2ā⋅m²⋅b̄ + 2ā⋅m⋅b̄ +2ē⋅ā⋅n - d̄⋅b̄⋅n) ⋅ (d̄ + 2ā⋅m + 2ā)
    c1 *= (m+1) ⋅ (d̄ + 3ā + 2ā⋅m) ⋅  (d̄ + ā + 2ā⋅m)

    c2 = -(m+2)
    c2a = -4ā²⋅c̄⋅m²
    c2a +=  ā⋅b̄²⋅m² + 2ā⋅b̄²⋅m - 4ā⋅c̄⋅m⋅d̄ - 8ā²⋅c̄⋅m + m⋅b̄²⋅d̄ - ā⋅ē^2 - d̄²⋅c̄ + b̄⋅ē⋅d̄ -4ā²⋅c̄
    c2a += -4ā⋅c̄⋅d̄ + ā⋅b̄² +b̄²⋅d̄
    c2 *= c2a ⋅ (ā⋅m + ā⋅n + ā + d̄) ⋅ (m + 1) ⋅ (d̄ + 2ā⋅m)

    return (c0, c1, c2)

end

## With an explicit connectioni
function connection(::Type{P}, q::Q) where
    {P <: Polynomials.AbstractPolynomial,
     Q<:Polynomials.AbstractPolynomial}

    n = degree(q)
    T = eltype(q)
    T = eltype(one(eltype(one(P)))*one(q)/1)
    cs = zeros(T, n+1)

    for k in 0:n
        for (i,val) in Connection{P,Q}(n, k)
            cs[1+k] = muladd(q[i], val, cs[1+k])
        end
    end
    X = Polynomials.indeterminate(P,q)
    ⟒(P){T,X}(cs)

end
