Rely on other packages for higher performance

Project.toml

@@ -4,19 +4,21 @@ authors = ["jverzani <jverzani@gmail.com>"]
 version = "0.1.0"
 
 [deps]
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
+HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
-
 [compat]
-julia = "1.0"
-Polynomials = "1.1"
-QuadGK = "1,2"
 Memoize = "0.4"
+Polynomials = "1.1.2"
+QuadGK = "1,2"
 SpecialFunctions = "0.09, 0.10"
+julia = "1.0"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -1,7 +1,7 @@
 # Special Polynomials
 
 [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://jverzani.github.io/SpecialPolynomials.jl/stable)
-![Docs from JuliaHub](https://juliahub.com/docs/SpecialPolynomials/)
+[![](https://img.shields.io/badge/docs-juliahub-blue.svg)](https://juliahub.com/docs/SpecialPolynomials/)
 [![Build Status](https://travis-ci.org/jverzani/SpecialPolynomials.jl.svg?branch=master)](https://travis-ci.org/jverzani/SpecialPolynomials.jl)
 
 A package providing various polynomial types (assuming different

docs/src/examples.md

@@ -81,7 +81,7 @@ julia> p3(u)
 ChebyshevU(- 0.125⋅U₁(x) + 0.3125⋅U₃(x))
 ```
 
-For most of the orthogonal polynomials, a conversion from the standard basis is provided, and a conversion between different parameter values for the same polynomial type are provded. Conversion methods between other polynomial types are not provided, but either evaluation, as above, or conversion through the `Polynomial` type is possible.
+For most of the orthogonal polynomials, a conversion from the standard basis is provided, and a conversion between different parameter values for the same polynomial type are provded. Conversion methods between other polynomial types are not provided, but either evaluation, as above, or conversion through the `Polynomial` type is possible. As possible, for the orthogonal polynomial types, conversion utilizes the `FastTransforms` package; this package can handle conversion between polynomials with very high degree.
 
 
 For the basis functions, the `basis` function can be used:
@@ -95,14 +95,22 @@ julia> h3(x)
 Polynomial(-12.0*x + 8.0*x^3)
 ```
 
-If the coefficients are known, they can be directly passed to the constructor:
+For numeric evaluation of just a basis polynomial of a classical orthogonal polynomial system, the `Basis` constructor provides a direct evaluation without the construction of an intermediate polynomial:
+
+```
+Basis(Hermite, 3)(0.5)
+```
+
+
+
+If the coefficients of a polynomial relative to the polynomial type are known, they can be directly passed to the constructor:
 
 ```jldoctest example
 julia> Laguerre{0}([1,2,3])
 Laguerre{0}(1⋅L₀(x) + 2⋅L₁(x) + 3⋅L₂(x))
 ```
 
-Some polynomial types are parameterized. The parameters are passed as in this example:
+Some polynomial types are parameterized, as above. The parameters are passed to the type, as in this example:
 
 ```jldoctest example
 julia> Jacobi{1/2, -1/2}([1,2,3])
@@ -144,6 +152,29 @@ julia> SpecialPolynomials.innerproduct(P, p4, p5)
 
 For each polynomial type, this package implements as many of the methods for polynomials defined in `Polynomials`, as possible. 
 
+### Evaluation
+
+Evalution, as seen, is done through making polynomial objects callable:
+
+```
+julia> P = Chebyshev
+Chebyshev
+
+julia> p = P([1,2,3,4])
+Chebyshev(1⋅T₀(x) + 2⋅T₁(x) + 3⋅T₂(x) + 4⋅T₃(x))
+
+julia> p(0.4)
+-4.016
+```
+
+By default, for classical orthogonal polynomials, the Clenshaw reduction formula is used. For such polynomials, an alternative is to use the hypergoemetric formulation. (The evaluation `Basis(P,n)(x)` uses this.) There is an unexported method to compute through this means:
+
+```
+julia> SpecialPolynomials.eval_hyper(P, coeffs(p), 0.4)
+-4.016
+```
+
+
 
 ### Arithemtic
 
@@ -227,6 +258,28 @@ julia> integrate(p, 0, 1)
 25//8
 ```
 
+### Conversion
+
+Expressing a polynomial in type `P` in type `Q` is done through several possible means:
+
+```jldoctest example
+julia> P,Q = Gegenbauer{1//3}, Gegenbauer{2//3}
+(Gegenbauer{1//3,T,N} where N where T, Gegenbauer{2//3,T,N} where N where T)
+
+julia> p = P([1,2,3.0])
+Gegenbauer{1//3}(1.0⋅Cᵅ₀(x) + 2.0⋅Cᵅ₁(x) + 3.0⋅Cᵅ₂(x))
+
+julia> convert(Q, p)
+Gegenbauer{2//3}(0.8⋅Cᵅ₀(x) + 1.0⋅Cᵅ₁(x) + 1.2000000000000002⋅Cᵅ₂(x))
+
+julia> p(variable(Q))
+Gegenbauer{2//3}(0.7999999999999999⋅Cᵅ₀(x) + 1.0⋅Cᵅ₁(x) + 1.1999999999999997⋅Cᵅ₂(x))
+
+julia> SpecialPolynomials._convert_cop(Q,p)
+Gegenbauer{2//3}(0.8⋅Cᵅ₀(x) + 1.0⋅Cᵅ₁(x) + 1.2⋅Cᵅ₂(x))
+```
+
+The first uses a method from the `FastTransforms` package. This package can handle polynomials of very high degree. It is used by default, as much as possible. The second uses polynomial evalution (Clenshaw evaluation) to perform the conversion. The third uses the structural equations for conversion, when possible, and defaults to converting through the `Polynomial` type
 
 ### Roots
 
@@ -312,15 +365,27 @@ julia> eigvals(SpecialPolynomials.jacobi_matrix(Legendre, 50 )) .|> isreal |> s
 (The roots of the classic orthogonal polynomials are all real and distinct.)
 
 
-The unexported `gauss_nodes_weights` function returns the nodes and weights. For many types (e.g., `Jacobie`, `Legendre`, `Hermite`, `Laguerre`) it uses an `O(n)` algorithm of Glaser, Liu, and Rokhlin, for others the `O(n²)` algorithm through the Jacobi matrix.
+The unexported `gauss_nodes_weights` function returns the nodes and weights. For many types (e.g., `Jacobie`, `Legendre`, `Hermite`, `Laguerre`). As possible, it uses the methods from the `FastGaussQuadratures` package, which provides `O(n)` algorithms, where the Jacobi matrix is `O(n²)`.
+
+
 
 ```jldoctest example
-julia> xs, ys = SpecialPolynomials.gauss_nodes_weights(Legendre{Float64}, 10)
-([-0.9739065285171717, -0.8650633666889845, -0.6794095682990244, -0.4333953941292472, -0.1488743389816312, 0.1488743389816312, 0.4333953941292472, 0.6794095682990244, 0.8650633666889845, 0.9739065285171717], [0.06667134430868804, 0.1494513491505808, 0.21908636251598188, 0.26926671930999635, 0.295524224714753, 0.295524224714753, 0.26926671930999635, 0.21908636251598188, 0.1494513491505808, 0.06667134430868804])
+julia> xs, ws = SpecialPolynomials.gauss_nodes_weights(Legendre, 4)
+([-0.8611363115940526, -0.3399810435848563, 0.3399810435848563, 0.8611363115940526], [0.34785484513745385, 0.6521451548625462, 0.6521451548625462, 0.34785484513745385])
+
+julia> basis(Legendre, 4).(xs)
+4-element Array{Float64,1}:
+ 1.1102230246251565e-16
+ -8.326672684688674e-17
+ -8.326672684688674e-17
+ 1.1102230246251565e-16
+
+julia> f(x) = x^7 - x^6; F(x) = x^8/8 - x^7/7;
+
+julia> sum(f(x)*w for (x,w) in zip(xs, ws)) - (F(1) - F(-1))
+5.551115123125783e-17
 ```
 
-!!! note
- For a broader and more robust implementation, the [FastGaussQuadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) package provides `O(n)` algorithms for many classic orthogonal polynomial types.
 
 ## Fitting
 
@@ -493,7 +558,6 @@ degree(p)
  The [ApproxFun](https://github.com/JuliaApproximation/ApproxFun.jl) package provides a framework to quickly and accuratately approximate functions using certain polynomial types. The choice of order and methods for most of Julia's built-in functions are conveniently provided.
 
 
-
 ## Plotting
 
 The `plot` recipe from the `Polynomials` package works as expected for

src/Orthogonal/Chebyshev.jl

@@ -40,14 +40,30 @@ basis_symbol(::Type{<:Chebyshev}) = "T"
 Polynomials.domain(::Type{<:Chebyshev}) = Polynomials.Interval(-1, 1, false, false)
 weight_function(::Type{<: Chebyshev}) = x -> one(x)/sqrt(one(x) - x^2)
 generating_function(::Type{<: Chebyshev}) = (t,x) -> (1-t*x)/(1-2*t*x - t^2)
+function classical_hypergeometric(P::Type{<:Chebyshev}, n, x) where {α}
+ as = (-n,n)
+ bs = (one(eltype(P))/2, )
+ pFq(as, bs, (1-x)/2)
+end
 
+function eval_basis(::Type{Chebyshev}, n, x)
+ if abs(x) <= 1
+ cos(n * acos(x))
+ elseif x > 1
+ cosh(n*acosh(x))
+ else
+ (-1)^n * cosh(n*acosh(-x))
+ end
+end
 
 abcde(::Type{<:Chebyshev}) = NamedTuple{(:a,:b,:c,:d,:e)}((-1, 0, 1, -1, 0))
 
+k0(P::Type{<:Chebyshev}) = one(eltype(P))
 k1k0(P::Type{<:Chebyshev}, n::Int) = iszero(n) ? one(eltype(P)) : 2*one(eltype(P))
 
-#kn(P::Type{<:Chebyshev}, n::Int) = iszero(n) ? one(eltype(P)) : (2*one(eltype(P)))^(n-1)
-#k1k_1(P::Type{<:Chebyshev}, n::Int) = n==1 ? 2*one(eltype(P)) : 4*one(eltype(P))
+norm2(P::Type{<:Chebyshev}, n::Int) = (iszero(n) ? 1 : 1/2) * pi
+ω₁₀(P::Type{<:Chebyshev}, n::Int) = iszero(n) ? one(eltype(P))/sqrt(2) : one(eltype(P)) # √(norm2(n+1)/norm2(n)
+
 
 leading_term(P::Type{<:Chebyshev}, n::Int) = iszero(n) ? one(eltype(P)) : (2*one(eltype(P)))^(n-1)
 
@@ -57,7 +73,7 @@ Bn(P::Type{<:Chebyshev}, n::Int) = zero(eltype(P))
 Cn(P::Type{<:Chebyshev}, n::Int) = one(eltype(P))
 
 # 
-C̃n(P::Type{<:Chebyshev}, ::Val{1}) = one(eltype(P))
+C̃n(P::Type{<:Chebyshev}, ::Val{1}) = one(eltype(P))/2
 ĉ̃n(P::Type{<:Chebyshev}, ::Val{0}) = one(eltype(P))/4
 ĉ̃n(P::Type{<:Chebyshev}, ::Val{1}) = Inf
 γ̃n(P::Type{<:Chebyshev}, n::Int) = (n==1) ? one(eltype(P))/2 : n*one(eltype(P))/4
@@ -182,12 +198,15 @@ function lagrange_barycentric_nodes_weights(::Type{<: Chebyshev}, n::Int)
  xs, ws
 end
 
-# noded for integrating against the weight function
-function gauss_nodes_weights(::Type{<:Chebyshev}, n::Int)
- xs = cos.(pi/2n * (2*(n:-1:1).-1))
- ws = pi/n * ones(n)
- xs, ws
-end
+# # noded for integrating against the weight function
+# function gauss_nodes_weights(::Type{<:Chebyshev}, n::Int)
+# xs = cos.(pi/2n * (2*(n:-1:1).-1))
+# ws = pi/n * ones(n)
+# xs, ws
+# end
+
+gauss_nodes_weights(p::Type{P}, n) where {P <: Chebyshev} =
+ FastGaussQuadrature.gausschebyshev(n)
 
 ##
 ## fitting coefficients
@@ -329,6 +348,15 @@ export MonicChebyshev
 @register_monic(MonicChebyshev)
 C̃n(P::Type{<:MonicChebyshev}, ::Val{1}) = one(eltype(P))
 
+
+##
+
+@register0 OrthonormalChebyshev AbstractCCOP0
+export OrthonormalChebyshev
+ϟ(::Type{<:OrthonormalChebyshev}) = Chebyshev
+ϟ(::Type{<:OrthonormalChebyshev{T}}) where {T} = Chebyshev{T}
+@register_orthonormal(OrthonormalChebyshev)
+
 ##
 ## --------------------------------------------------
 ##
@@ -424,6 +452,12 @@ ChebyshevU
 basis_symbol(::Type{<:ChebyshevU}) = "U"
 weight_function(::Type{<:ChebyshevU}) = x -> sqrt(one(x) - x^2)
 generating_function(::Type{<: ChebyshevU}) = (t, x) -> 1 / (1 - 2t*x + t^2)
+function classical_hypergeometric(P::Type{<:ChebyshevU}, n, x) where {α}
+ as = (-n,n+2)
+ bs = ((3one(eltype(P)))/2, )
+ (n+1)*pFq(as, bs, (1-x)/2)
+end
+
 Polynomials.domain(::Type{<:ChebyshevU}) = Polynomials.Interval(-1, 1)
 
 abcde(::Type{<:ChebyshevU}) = NamedTuple{(:a,:b,:c,:d,:e)}((-1, 0, 1, -3, 0))
@@ -432,6 +466,9 @@ kn(P::Type{<:ChebyshevU}, n::Int) = (2 * one(eltype(P)))^n
 k1k0(P::Type{<:ChebyshevU}, n::Int) = 2 * one(eltype(P))
 k1k_1(P::Type{<:ChebyshevU}, n::Int) = 4 * one(eltype(P))
 
+norm2(P::Type{<:ChebyshevU}, n::Int) = pi/2
+ω₁₀(P::Type{<:ChebyshevU}, n::Int) = one(eltype(P))
+
 
 # directly adding these gives a 5x speed up in polynomial evaluation
 An(::Type{<:ChebyshevU}, n::Int) = 2
@@ -512,3 +549,11 @@ export MonicChebyshevU
 @register_monic(MonicChebyshevU)
 
 
+@register0 OrthonormalChebyshevU AbstractCCOP0
+export OrthonormalChebyshevU
+ϟ(::Type{<:OrthonormalChebyshevU}) = ChebyshevU
+ϟ(::Type{<:OrthonormalChebyshevU{T}}) where {T} = ChebyshevU{T}
+@register_orthonormal(OrthonormalChebyshevU)
+
+
+

src/Orthogonal/Discrete/FallingFactorial.jl

@@ -78,6 +78,8 @@ function (p::FallingFactorial)(x)
  tot
 end
 
+k0(P::Type{<:FallingFactorial}) = one(eltype(P))
+
 Base.:*(p::FallingFactorial, q::FallingFactorial) = convert(FallingFactorial, convert(Polynomial,p)*convert(Polynomial,q))
 
 

src/Orthogonal/Gegenbauer.jl

@@ -29,12 +29,20 @@ Polynomials.domain(::Type{<:Gegenbauer{α}}) where {α} = Polynomials.Interval(-
 
 abcde(::Type{<:Gegenbauer{α}}) where {α} = NamedTuple{(:a,:b,:c,:d,:e)}((-1,0,1,-(2α+1),0))
 
-
+k0(P::Type{<:Gegenbauer}) = one(eltype(P))
 k1k0(P::Type{<:Gegenbauer{α}}, k::Int) where {α} =(one(eltype(P))*2*(α+k))/(k+1)
 
 function norm2(::Type{<:Gegenbauer{α}}, n) where{α}
  pi * 2^(1-2α) * gamma(n + 2α) / (gamma(n+1) * (n+α) * gamma(α)^2)
 end
+function ω₁₀(::Type{<:Gegenbauer{α}}, n) where{α}
+ val = gamma(n + 1 + 2α) / gamma(n + 2α)
+ val /= (n+1)
+ val *= (n + α)/(n+1+α)
+ sqrt(val)
+
+end
+
 weight_function(::Type{<:Gegenbauer{α}}) where {α} = x -> (1-x^2)^(α-1/2)
 generating_function(::Type{<:Gegenbauer{α}}) where {α} = (t,x) -> begin
  1/(1-2x*t +t^2)^α
@@ -57,3 +65,11 @@ b̂̃n(P::Type{<:Gegenbauer{α}}, n::Int) where {α} = zero(eltype(P))# one(S) *
 b̂̃n(P::Type{<:Gegenbauer{α}}, ::Val{N}) where {α,N} = zero(eltype(P))# one(S) * 4/α/(α+2) 
 ĉ̃n(P::Type{<:Gegenbauer{α}}, ::Val{0}) where {α} = zero(eltype(P)) #one(S) * 4/α/(α+2) 
 
+
+
+
+@registerN OrthonormalGegenbauer AbstractCCOP1 α
+export OrthonormalGegenbauer
+ϟ(::Type{<:OrthonormalGegenbauer{α}}) where {α} = Gegenbauer{α}
+ϟ(::Type{<:OrthonormalGegenbauer{α,T}}) where {α,T} = Gegenbauer{α,T}
+@register_orthonormal(OrthonormalGegenbauer)

src/Orthogonal/Hermite.jl

@@ -48,6 +48,7 @@ Polynomials.domain(::Type{<:Hermite}) = Polynomials.Interval(-Inf, Inf)
 
 abcde(::Type{<:Hermite}) = NamedTuple{(:a,:b,:c,:d,:e)}((1,0,0,-2,0))
 
+k0(P::Type{<:Hermite}) = one(eltype(P))
 function k1k0(P::Type{<:Hermite}, k::Int)
  val = 2*one(eltype(P))
  val
@@ -62,10 +63,8 @@ function classical_hypergeometric(::Type{<:Hermite}, n, x)
  (2x)^n * pFq(as, bs, -inv(x)^2)
 end
 
-function gauss_nodes_weights(P::Type{<:Hermite}, n)
- xs, ws = gauss_nodes_weights(ChebyshevHermite, n)
- xs/sqrt(2), ws/sqrt(2)
-end
+gauss_nodes_weights(p::Type{P}, n) where {P <: Hermite} =
+ FastGaussQuadrature.gausshermite(n)
 
 
 ## Overrides