diff --git a/Project.toml b/Project.toml
index 4a38135..3843485 100644
--- a/Project.toml
+++ b/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"
diff --git a/README.md b/README.md
index f35eb6f..1d6044c 100644
--- a/README.md
+++ b/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
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 339099b..6961abe 100644
--- a/docs/src/examples.md
+++ b/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
diff --git a/src/Orthogonal/Chebyshev.jl b/src/Orthogonal/Chebyshev.jl
index d19f109..1a4b166 100644
--- a/src/Orthogonal/Chebyshev.jl
+++ b/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)
+
+
+
diff --git a/src/Orthogonal/Discrete/FallingFactorial.jl b/src/Orthogonal/Discrete/FallingFactorial.jl
index bebe482..37e5854 100644
--- a/src/Orthogonal/Discrete/FallingFactorial.jl
+++ b/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))
 
 
diff --git a/src/Orthogonal/Gegenbauer.jl b/src/Orthogonal/Gegenbauer.jl
index 61d2512..34130b7 100644
--- a/src/Orthogonal/Gegenbauer.jl
+++ b/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)
diff --git a/src/Orthogonal/Hermite.jl b/src/Orthogonal/Hermite.jl
index 601c86d..ebdac53 100644
--- a/src/Orthogonal/Hermite.jl
+++ b/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
diff --git a/src/Orthogonal/Jacobi.jl b/src/Orthogonal/Jacobi.jl
index a0ce6d4..0cc1710 100644
--- a/src/Orthogonal/Jacobi.jl
+++ b/src/Orthogonal/Jacobi.jl
@@ -31,7 +31,7 @@ basis_symbol(::Type{<:Jacobi{α,β}}) where {α,β} = "Jᵅᵝ"
 Polynomials.domain(::Type{<:Jacobi{α, β}}) where {α, β} = Polynomials.Interval(-1, 1, β >= 0, α >= 0)
 abcde(::Type{<:Jacobi{α,β}})  where {α,β} = NamedTuple{(:a,:b,:c,:d,:e)}((-1,0,1,-(α+β+2),β-α))
 
-
+k0(P::Type{<:Jacobi}) = one(eltype(P))
 function k1k0(P::Type{<:Jacobi{α,β}}, n::Int) where {α,β}
     n == -1  &&  return one(eltype(P))/1 
     γ = 2n + α + β
@@ -68,20 +68,8 @@ generating_function(::Type{<:Jacobi{α, β}}) where {α, β} = (t,x) -> begin
     2^(α + β) * 1/R  * (1 - t + R)^(-α) * (1 + t + R)^(-β)
 end
 
-function  gauss_nodes_weights(P::Type{<:Jacobi{α,β}}, n)  where {α,β}
-    # we don't  have a  good starting point  unless α=β
-    if α == β
-        xs, ws = glaser_liu_rokhlin_gauss_nodes(basis(MonicJacobi{α,β},n))
-        λ = kn(P,n)^2
-        xs, ws/λ
-    else
-        J = jacobi_matrix(P, n)
-        eig = eigen(J, extrema(P)...)
-        wts =  π̃βn(P,0) * (eig.vectors[1,:]).^2
-        eig.values,  wts
-    end
-
-end
+gauss_nodes_weights(p::Type{P}, n) where {α, β, P <: Jacobi{α, β}} =
+    FastGaussQuadrature.gaussjacobi(n, α, β)
 
 function classical_hypergeometric(::Type{<:Jacobi{α, β}}, n, x) where {α,β}
 
@@ -99,6 +87,14 @@ function norm2(::Type{<:Jacobi{α, β}}, n) where{α, β}
     α > -1 && β > -1 || throw(ArgumentError("α, β > -1 is necessary"))
     2^(α+β+1)/(2n+α+β+1) * (Γ(n+α+1) *  Γ(n + β +1))/(Γ(n+α+β+1)*Γ(n+1))
 end
+function ω₁₀(::Type{<:Jacobi{α, β}}, n) where{α, β}
+    val = Γ(n+1+α+1)/Γ(n+α+1)
+    val *= Γ(n + 1+ β +1) / Γ(n + β +1)
+    val *= Γ(n+α+β+1) / Γ(n+1+α+β+1)
+    val *= Γ(n+1) / Γ(n+1+1)
+    val *= (2n+α+β+1) / (2(n+1)+α+β+1)
+    sqrt(val)
+end
 
 # overrides
 B̃n(P::Type{<:Jacobi{α,β}}, ::Val{0}) where {α, β} = iszero(α+β) ? (α-β)*one(eltype(P))/2 : (α-β)*one(eltype(P))/(2(α+β+2))
@@ -122,18 +118,12 @@ export MonicJacobi
 
 @register_monic(MonicJacobi)
 
-pqr_start(P::Type{MonicJacobi{α,α}}, n) where {α} =  0.0 # zero(eltype(P))
-pqr_symmetry(P::Type{<:MonicJacobi{α,α}}) where {α} = true
-function pqr_weight(P::Type{<:MonicJacobi{α,α}}, n, x, dπx) where {α}
-    # https://www.chebfun.org/publications/HaleTownsend2013a.pdf
-    # Eqn  (1.4)
-    # Cnαβ should be  computed  using asymptotic formula  for larger n (§3.2.3)
-    # XXX TThere i ss ome constant  that  makes this not work....
-    β = α
-    Cnαβ = 2^(α+β+1)  
-    Cnαβ *= gamma(n + α +  1) * gamma(n + β + 1) / gamma(n + α + β + 1)
-    Cnαβ /= gamma(1 + n)
-    val = Cnαβ / (1-x^2) / dπx^2
-    val
-end
-                      
+
+@registerN OrthonormalJacobi AbstractCCOP2 α β
+export OrthonormalJacobi
+ϟ(::Type{<:OrthonormalJacobi{α,β}}) where {α,β} = Jacobi{α,β}
+ϟ(::Type{<:OrthonormalJacobi{α,β,T}}) where {α,β,T} = Jacobi{α,β,T}
+
+@register_orthonormal(OrthonormalJacobi)
+
+
diff --git a/src/Orthogonal/Laguerre.jl b/src/Orthogonal/Laguerre.jl
index 7101951..4b845d0 100644
--- a/src/Orthogonal/Laguerre.jl
+++ b/src/Orthogonal/Laguerre.jl
@@ -50,20 +50,23 @@ Polynomials.domain(::Type{<:Laguerre}) = Polynomials.Interval(0, Inf)
 
 abcde(::Type{<:Laguerre{α}})  where {α} = NamedTuple{(:a,:b,:c,:d,:e)}((0,1,0,-1,α+1))
 
-k1k0(P::Type{<:Laguerre{α}}, k::Int) where {α} = k <= 0  ? -one(eltype(P)) : -one(eltype(P))/(k+1)
+k0(P::Type{<:Laguerre}) = one(eltype(P))
+k1k0(P::Type{<:Laguerre{α}}, k::Int) where {α} = -one(eltype(P))/(k+1) # k >=0
 
 
-# function kn(P::Type{<:Laguerre{α}}, n::Int)  where {α}
-#     (-one(eltype(P)))^n/factorial(n) #Γ(1+n)
-# end
-
-#k1k_1(P::Type{<:Laguerre{α}}, k::Int) where {α} =  k <= 0  ? zero(eltype(P)) : one(eltype(P))/(k+1)/k
-
 function norm2(::Type{<:Laguerre{α}}, n) where{α}
     iszero(α) && return one(α)
-    gamma(n + α + 1) / gamma(n+1)
+    Γ(n + α + 1) / Γ(n+1)
+end
+
+function ω₁₀(::Type{<:Laguerre{α}}, n) where{α}
+    iszero(α) && return 1.0
+    val = Γ(n + 1 + α + 1) / Γ(n + α + 1)
+    val /= n + 1 # Γ(n+1)/Γ(n+2)
+    sqrt(val)
 end
 
+
 weight_function(::Type{<:Laguerre{α}}) where {α} = x  -> x^α * exp(-x)
 generating_function(::Type{<:Laguerre{α}}) where {α} = (t,x) -> begin
     exp(-t*x/(1-t))/(1-t)^(α-1)
@@ -75,12 +78,8 @@ function classical_hypergeometric(::Type{<:Laguerre{α}}, n, x) where {α}
     Pochhammer_factorial(α+1,n)*pFq(as, bs, x)
 end
 
-function gauss_nodes_weights(P::Type{<:Laguerre{α}}, n) where {α}
-    α <= -1 && throw(ArgumentError("α  too small"))
-    xs,  ws  =  glaser_liu_rokhlin_gauss_nodes(basis(MonicLaguerre{α},n))
-    λ = kn(P,n)^2
-    xs, ws/λ
-end
+gauss_nodes_weights(p::Type{P}, n) where {α, P <: Laguerre{α}} =
+    FastGaussQuadrature.gausslaguerre(n,α)
 
         
 ## Overrides
@@ -119,12 +118,9 @@ export MonicLaguerre
 ϟ(::Type{<:MonicLaguerre{α,T}}) where {α,T} = Laguerre{α,T}
 @register_monic(MonicLaguerre)
 
-#
-pqr_start(P::Type{MonicLaguerre{α}}, n) where {α} = 2/(4n+2α+2)
-pqr_symmetry(P::Type{<:MonicLaguerre{α}})  where {α} = false
-function pqr_weight(P::Type{<:MonicLaguerre{α}}, n, x, dπx) where {α}
-    val =  one(eltype(P)) 
-    val /= (x*dπx^2)
-    val = 1/(x*dπx^2)    
-    val
-end
+
+@registerN OrthonormalLaguerre AbstractCCOP1 α
+export OrthonormalLaguerre
+ϟ(::Type{<:OrthonormalLaguerre{α}}) where {α} = Laguerre{α}
+ϟ(::Type{<:OrthonormalLaguerre{α,T}}) where {α,T} = Laguerre{α,T}
+@register_orthonormal(OrthonormalLaguerre)
diff --git a/src/Orthogonal/Legendre.jl b/src/Orthogonal/Legendre.jl
index 8c0d027..9fc5c3d 100644
--- a/src/Orthogonal/Legendre.jl
+++ b/src/Orthogonal/Legendre.jl
@@ -50,21 +50,21 @@ basis_symbol(::Type{<:Legendre})  = "P"
 Polynomials.domain(::Type{<:Legendre}) = Polynomials.Interval(-1, 1)
 
 
-#kn(P::Type{<:Legendre}, n::Int)  = kn(Gegenbauer{1/2, eltype(P)}, n)
+k0(P::Type{<:Legendre}) = one(eltype(P))
 k1k0(P::Type{<:Legendre}, n::Int) = (one(eltype(P))*(2n+1))/(n+1) #k1k0(Gegenbauer{1/2, eltype(P)}, n)
-#k1k_1(P::Type{<:Legendre}, n::Int) = k1k_1(Gegenbauer{1/2, eltype(P)}, n)
 
 norm2(::Type{<:Legendre}, n) = 2/(2n+1)
+ω₁₀(P::Type{<:Legendre}, n) = sqrt((one(eltype(P))*(2n+1))/(2n+3))
 weight_function(::Type{<:Legendre})  = x -> one(x)
 generating_function(::Type{<:Legendre}) = (t, x)  -> 1/sqrt(1 - 2x*t +t^2)
-
-# gauss nodes
-function gauss_nodes_weights(P::Type{<:Legendre}, n)
-    xs,  ws =  glaser_liu_rokhlin_gauss_nodes(basis(MonicLegendre,n))
-    λ = kn(P,n)^2
-    xs, ws/λ
+function classical_hypergeometric(::Type{<:Legendre}, n, x)
+    as = (-n, n+1)
+    bs = (1, )
+    pFq(as, bs, (1-x)/2)
 end
 
+gauss_nodes_weights(p::Type{P}, n) where {P <: Legendre} =
+    FastGaussQuadrature.gausslegendre(n)
 
 # overrides
 
@@ -145,3 +145,9 @@ export MonicShiftedLegendre
 @register_monic(MonicShiftedLegendre)
 
 
+@register0 OrthonormalLegendre AbstractCCOP0
+export OrthonormalLegendre
+ϟ(::Type{<:OrthonormalLegendre}) = Legendre
+ϟ(::Type{<:OrthonormalLegendre{T}}) where {T} = Legendre{T}
+@register_orthonormal(OrthonormalLegendre)
+
diff --git a/src/Orthogonal/abstract.jl b/src/Orthogonal/abstract.jl
index bfaa702..be1dee7 100644
--- a/src/Orthogonal/abstract.jl
+++ b/src/Orthogonal/abstract.jl
@@ -101,7 +101,7 @@ macro registerN(name,  parent, params...)
             $poly{$(αs...),promote_type(T,S)}
 
         $poly{$(αs...),T}(n::Number, var::Polynomials.SymbolLike = :x) where {$(αs...),T} =
-            $poly{$(αs...)}(T[n], var)
+            n * one($poly{$(αs...),T}, var)
         $poly{$(αs...)}(n::S, var::Polynomials.SymbolLike = :x)        where {$(αs...), S<:Number} =
             $poly{$(αs...),S}(n,var)
         $poly{$(αs...),T}(var::Polynomials.SymbolLike=:x)              where {$(αs...), T} =
@@ -147,6 +147,31 @@ macro register_monic(name)
     end
 end
 
+## Make ortho
+macro register_orthonormal(name)
+
+    orthonormal=esc(name)
+    
+    quote
+        SpecialPolynomials.isorthonormal(::Type{<:$orthonormal}) = true
+        SpecialPolynomials.abcde(P::Type{<:$orthonormal})  = abcde(ϟ(P))
+        SpecialPolynomials.basis_symbol(P::Type{<:$orthonormal})  = basis_symbol(ϟ(P)) * "̃"
+        SpecialPolynomials.weight_function(P::Type{<:$orthonormal}) = weight_function(ϟ(P))
+        Polynomials.domain(P::Type{<:$orthonormal}) = domain(ϟ(P))
+        SpecialPolynomials.k0(::Type{P}) where{P <: $orthonormal} = k0(ϟ(P)) / sqrt(norm2(ϟ(P),0)) 
+        SpecialPolynomials.k1k0(::Type{P},  n::Int) where{P <: $orthonormal} = k1k0(ϟ(P),n) /  ω₁₀(ϟ(P), n)
+        SpecialPolynomials.ω₁₀(P::Type{<:$orthonormal}, n::Int) = one(eltype(P))
+        
+        SpecialPolynomials.B̃n(P::Type{<:$orthonormal}, n::Int) =  B̃n(ϟ(P),n)
+        SpecialPolynomials.B̃n(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  B̃n(ϟ(P),v)
+        SpecialPolynomials.C̃n(P::Type{<:$orthonormal}, n::Int) =  C̃n(ϟ(P),n)
+        SpecialPolynomials.C̃n(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  C̃n(ϟ(P),v)
+        SpecialPolynomials.ẫn(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  ẫn(ϟ(P),v)
+        SpecialPolynomials.b̂̃n(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  b̂̃n(ϟ(P),v)
+        SpecialPolynomials.ĉ̃n(P::Type{<:$orthonormal}, v::Val{N}) where {N} =  ĉ̃n(ϟ(P),v)
+    end
+end
+
 # P̃ = P(αx+β), so P̃'' = α^2P''(α x+ β); P̃' = αP'(αx+β)
 macro register_shifted(name, alpha, beta)
     shifted = esc(name)
diff --git a/src/Orthogonal/ccop.jl b/src/Orthogonal/ccop.jl
index c681f3c..67cbcc0 100644
--- a/src/Orthogonal/ccop.jl
+++ b/src/Orthogonal/ccop.jl
@@ -161,13 +161,13 @@ end
 
 # An, Bn,  Cn
 # p_{n+1} = (An*x + Bn)⋅p_n + Cn⋅p_{n-1}
-An(P::Type{<:AbstractCOP}, n::Int) = Ãn(P,n) * k1k0(P, n)
+An(P::Type{<:AbstractCOP}, n::Int) = Ãn(P,n) * k1k0(P, n) 
 function  Ãn(P::Type{<:AbstractCOP}, n::Int)
     a,b,c,d,e = abcde(P)
     Ãn(P, a,b,c,d,e ,n) 
 end
 
-Bn(P::Type{<:AbstractCOP}, n::Int) = B̃n(P,n) * k1k0(P,n)
+Bn(P::Type{<:AbstractCOP}, n::Int) = B̃n(P,n) * k1k0(P,n) 
 function B̃n(P::Type{<:AbstractCOP}, n::Int)
     a,b,c,d,e = abcde(P)
     B̃n(P, a,b,c,d,e ,n) 
@@ -366,13 +366,26 @@ end
 ##
 ## Conversion
 
+
 function Base.convert(::Type{Q},  p::P)  where {Q <: Polynomials.StandardBasisPolynomial, P <: AbstractCOP} 
     p(variable(Q, p.var))
 end
 function Base.convert(::Type{Q},  p::P)  where {Q <: AbstractCCOP,  P <: Polynomials.StandardBasisPolynomial}
     _convert_cop(Q, p)
 end
+
+## Conversion
+## * use FastTransforms, when available for T <: AbstractFloat; see connection.jl
+## * use  _convert_cop when possible (needs to match σ)
+## * use  convesion  through Polynomial type
 function Base.convert(::Type{Q}, p::P)  where  {Q <: AbstractCCOP,  P <: AbstractCCOP}
+    _convert(Q, p)
+end
+
+# work around method ambiguity introducted in abstract
+# dispatch  to  specific FastTransform  method  (defined in `connection.jl`) or
+# use this default
+function  _convert(::Type{Q}, p::P) where {Q <: AbstractCCOP,  P <: AbstractCCOP}
 
     a,b,c,d,e = abcde(P)
     ā,b̄,c̄,d̄,ē = abcde(Q)
@@ -383,45 +396,23 @@ function Base.convert(::Type{Q}, p::P)  where  {Q <: AbstractCCOP,  P <: Abstrac
         T = eltype(Q)
         convert(Q, convert(Polynomial{T}, p))
     end
-end
+end    
 
 ##
 ## --------------------------------------------------
 # # scalar ops
-# function Base.:+(p::P, c::S) where {P <: AbstractCOP, S<:Number}
-#     R = promote_type(eltype(p), S)
-#     as = R[a for a in coeffs(p)]
-#     if length(as) >= 1
-#         as[1] += c  # *assume* basis(P,0) = 1
-#     else
-#         push!(as, c)
-#     end
-#     ⟒(P)(as, p.var)
-# end
 
 #  avoid dispatch when N is known
 function Base.:+(p::P, c::S) where {T, N, P<:AbstractCOP{T,N}, S<:Number}
-    R = promote_type(T,S)
-    iszero(c) && return (N == 0 ? zero(⟒(P){R},p.var) :  ⟒(P){R,N}(R.(p.coeffs), p.var))
-    N == 0 && return ⟒(P)(R[c], p.var)
-    N == 1 && return ⟒(P)(R[p[0]+c], p.var)
-    cs = R[iszero(i) ? p[i]+c : p[i] for i in 0:N-1]
+    c′ =  c / k0(P) #one(T) * ⟒(P)(c)[0] #  c / k0(P), needed to add to a coefficient
+    R = promote_type(promote_type(T,S), typeof(inv(k0(P))))
+    iszero(c) && return (N == 0 ? zero(⟒(P){R}, p.var) :  ⟒(P){R,N}(R.(p.coeffs), p.var))
+    N == 0 && return ⟒(P)(R(c), p.var)
+    N == 1 && return ⟒(P)(R[p[0] + c′], p.var)
+    cs = R[iszero(i) ? p[0]+c′ : p[i] for i in 0:N-1]
     return ⟒(P){R,N}(cs, p.var)
 end
 
-# function Base.:*(p::P, c::S) where {T, N, P<:AbstractCOP{T,N}, S <: Number}
-#     R = promote_type(T,S)
-#     iszero(c)  && return zero(⟒(P){R})
-#     ⟒(P){R,N}(p.coeffs .* c, p.var)
-# end
-
-# function Base.:/(p::P, c::S) where {T,N,P<:AbstractCOP{T,N},S <: Number}
-#     R = eltype(one(T)/one(S))
-#     isinf(c)  && return zero(⟒(P){R})
-#     ⟒(P){R,N}(p.coeffs ./ c, p.var)
-# end
-
-
 ##
 ## multiply/addition/divrem with P{α…,T,N} =  Q{α...,T, M}
 ## We don't have (p::P{N},q::P{M}) where  {N,M, P<:AbstractCOP} as an available signature
@@ -430,9 +421,10 @@ function ⊕(p::P, q::Q) where {T,N,S,M, P <: AbstractCOP{T,N}, Q <: AbstractCOP
 
     #@assert  ⟒(P) == ⟒(Q)
     #@assert eltype(p) == eltype(q)
+    R = promote_type(T,S)
+    Polynomials.isconstant(p)  && return q + p(0) 
+    Polynomials.isconstant(q)  && return p + q(0) 
 
-    Polynomials.isconstant(p)  && return q + p[0]
-    Polynomials.isconstant(q)  && return p + q[0]    
     p.var != q.var && throw(ArgumentError("Variables don't  match"))    
 
     if N==M
@@ -468,12 +460,14 @@ end
 
 function ⊗(p::P, q::Q) where {P <: AbstractCOP, Q <: AbstractCOP}
 
-    Polynomials.isconstant(p)  && return q * p[0]
-    Polynomials.isconstant(q)  && return p * q[0]    
+    Polynomials.isconstant(p)  && return q * p(0)
+    Polynomials.isconstant(q)  && return p * q(0)
     p.var != q.var && throw(ArgumentError("Variables don't  match"))    
 
     # use connection for linearization;  note:  evalauation  is  faster than _convert_cop
-    convert(⟒(P), convert(Polynomial, p) * convert(Polynomial, q))
+    p′,q′ = _convert_cop.(Polynomial, (p,q))
+    _convert_cop(⟒(P), p′ * q′)
+#    convert(⟒(P), convert(Polynomial, p) * convert(Polynomial, q))
 
 end
 
@@ -536,7 +530,7 @@ function Polynomials.derivative(p::P, order::Integer=1) where {P <:AbstractCOP}
             ps[1+n-1] -= pn*c/a
         end
     end
-    a,b,c = ân(P,0),b̂n(P,0),ĉn(P,0)
+    a,b = ân(P,0),b̂n(P,0)
     p1 = ps[1+1]
     as[1+0] = p1/a
 
@@ -550,34 +544,37 @@ function Polynomials.integrate(p::P, C::Number=0) where {P <: AbstractCOP}
     T,S = eltype(p), typeof(C)
     R = promote_type(typeof(one(T) / 1), S)
     Q = ⟒(P){R}
-    #if hasnan(p) || hasnan(C)
-    #    error("XXX nan")
-    #end
+    if hasnan(p) || isnan(C)
+        return   Q(NaN)
+    end
+
     n = degree(p)
-    if n == 0
-        return Q([C, p[0]], p.var)
+    if n == -1
+        return zero(Q, p.var)
+    elseif n == 0
+        return C*one(Q, p.var)  + p(0)*variable(Q, p.var)
     end
     
     as = zeros(R, n + 2)
 
     # case d=0 we do by hand,  as
-    # P_0(x) = 1, so ∫P_o = x = variable(P)
+    # P_0(x) = c, so ∫P_o = x = c*variable(P)
     c₀,c₁ = coeffs(variable(p))
     pd = first(p.coeffs)
-    as[1] = pd*c₀
-    as[2] = pd*c₁
+    as[1] = pd * c₀
+    as[2] = pd * (c₁ * k0(Q))
     @inbounds for d in 1:n
         pd = p.coeffs[d+1]
         as[1 + d + 1] += pd * ân(Q, d)
         as[1 + d]     += pd * b̂n(Q, d)
-        if  d > 0
+        if  d > 1
             as[1 + d - 1] += pd * ĉn(Q, d)
         end
     end
 
     # adjust constant
     ∫p = Q(as,  p.var)
-    return ∫p + (R(C) - ∫p(0))
+    return ∫p - ∫p(0) + Q(C) 
 
 end
 
diff --git a/src/Orthogonal/connection.jl b/src/Orthogonal/connection.jl
index 025ed89..2d5fb78 100644
--- a/src/Orthogonal/connection.jl
+++ b/src/Orthogonal/connection.jl
@@ -48,23 +48,23 @@ function _convert_cop(::Type{Q}, p::P) where {P <: ConvertibleTypes,
                                                Q <: ConvertibleTypes}
 
     d = degree(p)
-    T,S = eltype(one(Q)), eltype(p)  #
+    T,S = eltype(Q), eltype(p)  #
     R = typeof(one(promote_type(T,S))/1)
 
     as = zeros(R, 1+d)
 
-    λⱼⱼ = one(R) #  kn(P,0,R)/kn(Q,0,R) = 1
+    λⱼⱼ = one(R)  * k0(P) / k0(Q)
     
     for j in  0:d
 
-        λ = λⱼⱼ
+        λ = λⱼⱼ 
 
-        λⱼⱼ *= k1k0(P,j) / k1k0(Q,j)
+        λⱼⱼ *= k1k0(P,j) / k1k0(Q,j)  
 
         pⱼ = p[j]
         iszero(pⱼ) && continue
 
-        C̃ʲⱼ = one(R)
+        C̃ʲⱼ = one(R)  
 
         as[1+j] += pⱼ * (λ * C̃ʲⱼ)
 
@@ -85,6 +85,135 @@ function _convert_cop(::Type{Q}, p::P) where {P <: ConvertibleTypes,
     ⟒(Q){R}(as, p.var)
 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 {
diff --git a/src/Orthogonal/cop.jl b/src/Orthogonal/cop.jl
index 469362b..1050959 100644
--- a/src/Orthogonal/cop.jl
+++ b/src/Orthogonal/cop.jl
@@ -45,23 +45,29 @@ abcde(::Type{<:AbstractCOP}) = throw(ArgumentError("No default method"))
 """
     k1k0
 
-Let  `kᵢ` be the leading  coeffiecient of the  polynomial  in  the standard basis. This  function implement  `k₍ᵢ₊₁₎/kᵢ` for `i ≥  0`.
+Let  `kᵢ` be the leading  coeffiecient of the  polynomial  in  the standard basis. 
+This  function implements  `k₍ᵢ₊₁₎/kᵢ` for `i ≥ 0`.
 
-With  the assumption that  `k₀ = 1`, the values `kᵢ` and `k₍ᵢ₊₁₎/k₍ᵢ₋₁₎` can  be generated.
+The values `kᵢ` and `k₍ᵢ₊₁₎/k₍ᵢ₋₁₎` can  be generated.
 
 The default value  leads to monic polynomials.
 """
 k1k0(::Type{P},  i::Int) where {P<:AbstractCOP} =  one(eltype(P))
+k1k0(::Type{P},  i::Int) where {P<:Polynomials.StandardBasisPolynomial} =  one(eltype(P))
+
+"""
+      k0(::Type{P})
+
+This is  `basis(P,0)`,  most  often  just `1`, but may be different with some normalizations, such as `Orthonormal`.
+"""
+k0(::Type{P})  where{P <: AbstractCOP}  =  one(eltype(P))
+k0(::Type{P}) where {P <: Polynomials.StandardBasisPolynomial} = one(eltype(P))
 
 # Set defaults to be monic
 # For non-monic subtypes of `AbstractCOP` only need k1k0  to be defined, as `kn` and `k1k_1` come for free
 
-# leading term (Section 3 table)  is kn
-leading_term(P::Type{<:AbstractCOP},  n::Int) =  kn(P, n)
-
-# kn = prod(k1k0(i) for i in 0:n-1)  *  kn(P,0)
-# **Assume** basis(P,0) = 1, so kn(P,0)  = 1
-kn(::Type{P},  n::Int) where{P <: AbstractCOP} = iszero(n) ?  one(eltype(P)) : prod(k1k0(P,i) for i in  0:n-1)
+# kn = prod(k1k0(i) for i in 0:n-1)  *  k0(P)
+kn(::Type{P},  n::Int) where{P <: AbstractCOP} = foldr(*, (k1k0(P,i) for i in 0:n-1), init=k0(P)) 
 
 # k₍ᵢ₊₁₎/k₍ᵢ₋₁₎ =  (k₍ᵢ₊₁₎/kᵢ) ⋅ (kᵢ/k₍ᵢ₋₁₎)
 function k1k_1(::Type{P},  i::Int) where {P<:AbstractCOP}
@@ -69,6 +75,13 @@ function k1k_1(::Type{P},  i::Int) where {P<:AbstractCOP}
     k1k0(P,i)*k1k0(P,i-1)
 end
 
+# leading term (Section 3 table)  is kn
+leading_term(P::Type{<:AbstractCOP},  n::Int) =  kn(P, n)
+
+# square root of ratio of norm2(P,n+1)/norm2(P,n)
+# Let ωᵢ = √{∫ πᵢ² dw}, this is ωᵢ₊₁ ÷ ωᵢ
+ω₁₀(::Type{P},n)  where {P <:  AbstractCOP} = sqrt(norm2(P,n+1)/norm2(P,n))
+
 
 # Can't  change  N  here
 function Base.setindex!(p::AbstractCOP{T,N}, value::Number, idx::Int) where {T, N}
@@ -87,26 +100,50 @@ end
 Polynomials.degree(p::AbstractCOP{T,N})  where {T,N} = N-1
 Polynomials.isconstant(p::AbstractCOP) = degree(p) <=  0
 
+##  Evaluation
 
 """
     Clenshaw evaluation of an orthogonal polynomial 
 """
 function eval_cop(P::Type{<:AbstractCOP{T,N}}, cs, x::S) where {T,N,S}
+    N == 0 && return zero(T) * zero(S)
+    N == 1 && return (cs[1] * k0(P)) *  one(S)
+    _eval_cop(P,cs, x)
+end
+
+function _eval_cop(P::Type{<:AbstractCOP{T,N}}, cs, x::S) where {T,N,S}
     if @generated
-        N == 0 && return zero(T) * zero(S)
-        N == 1 && return cs[1] * one(S)
-        #SS = eltype(one(S))
-        Δ0 = :(cs[N-1])
-        Δ1 = :(cs[N])
-        for i in N-1:-1:2
-            a = :(cs[i - 1] - Δ1 * Cn(P, i-1))
-            b = :(Δ0 + Δ1 * muladd(x, An(P,i-1),Bn(P,i-1)))
-            Δ0 = :(a)
-            Δ1 = :(b)
+        quote
+            Δ0 = cs[end - 1]
+            Δ1 = cs[end]
+            @inbounds for i in N-1:-1:2
+                Δ0, Δ1 = cs[i - 1] - Δ1 * Cn(P, i-1), Δ0 + Δ1 * muladd(x, An(P,i-1), Bn(P,i-1))
+            end
+            p₀ = k0(P)
+            p₁ =  muladd(x, An(P,0),  Bn(P,0)) * p₀
+            Δ0 * p₀  + Δ1 * p₁
         end
-        Δ0 + Δ1* muladd(x, An(P,0), Bn(P,0))
     else
         clenshaw_eval(P, cs, x)
     end
 end
 
+#  evaluate  basis vector through hypergeometric formulation
+(B::Basis{P})(x) where  {P <: AbstractCOP} = eval_basis(P, B.n, x)
+
+# Evaluate a basis vector without realizing it and without using Clenshaw
+eval_basis(::Type{P}, n, x) where {P <: AbstractCOP} = classical_hypergeometric(P, n, x)
+
+"""
+    eval_hyper(::P, cs,  x::S)
+
+Evaluate polynomial `P(cs)` by  computing value for each basis vector from its hypergeometric representation
+"""
+function eval_hyper(P::Type{<:AbstractCOP}, cs, x::S) where {S}
+    isempty(cs) &&  return zero(S)
+    tot = cs[1] * one(S)
+    for i in 2:length(cs)
+        tot += cs[i] * Basis(P, i-1)(x)
+    end
+    tot
+end
diff --git a/src/Orthogonal/glaser-liu-rokhlin.jl b/src/Orthogonal/glaser-liu-rokhlin.jl
index d4e6b72..3ca0904 100644
--- a/src/Orthogonal/glaser-liu-rokhlin.jl
+++ b/src/Orthogonal/glaser-liu-rokhlin.jl
@@ -5,8 +5,43 @@
 
 Find Gauss nodes and weights of the  basis  polynomial  `π = basis(P,n)`. The nodes  are the  roots of the  polynomial and the weights are computed  from known formulas.
 
-The method from  [Glaser,  Liu, and Rokhlin](DOI: 10.1137/06067016X) is  used. This is  an  O(n)  method (whereas the method  basedon the Jacobi  matrix is O(n^2)). This method fits easily  into the framework  provided  through  the `AbstractCCOP` types.  The [FastGaussQuadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) package provides even more efficient  algorithms and  pays attention to numerical  issues, examples of  which  can  be found in  [Hale and Townsend](https://core.ac.uk/download/pdf/9403469.pdf). The `FastGuassQuadrature` package is recommended for actual  use  of these values. 
+The method from  [Glaser,  Liu, and Rokhlin](DOI: 10.1137/06067016X) is  used. This is  an  O(n)  method (whereas the method  basedon the Jacobi  matrix is O(n^2)). This method fits easily  into the framework  provided  through  the `AbstractCCOP` types.  The [FastGaussQuadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) package provides even more efficient  algorithms and  pays attention to numerical  issues, examples of  which  can  be found in  [Hale and Townsend](https://core.ac.uk/download/pdf/9403469.pdf). The `FastGuassQuadrature` package is used when  available.
+
+
+An example:
+
+```
+function  gauss_nodes_weights(P::Type{<:Jacobi{α,β}}, n)  where {α,β}
+    # we don't  have a  good starting point  unless α=β
+    if α == β
+        xs, ws = glaser_liu_rokhlin_gauss_nodes(basis(MonicJacobi{α,β},n))
+        λ = kn(P,n)^2
+        xs, ws/λ
+    else
+        J = jacobi_matrix(P, n)
+        eig = eigen(J, extrema(P)...)
+        wts =  π̃βn(P,0) * (eig.vectors[1,:]).^2
+        eig.values,  wts
+    end
 
+end
+    
+
+pqr_start(P::Type{MonicJacobi{α,α}}, n) where {α} =  0.0 # zero(eltype(P))
+pqr_symmetry(P::Type{<:MonicJacobi{α,α}}) where {α} = true
+function pqr_weight(P::Type{<:MonicJacobi{α,α}}, n, x, dπx) where {α}
+    # https://www.chebfun.org/publications/HaleTownsend2013a.pdf
+    # Eqn  (1.4)
+    # Cnαβ should be  computed  using asymptotic formula  for larger n (§3.2.3)
+    # XXX TThere i ss ome constant  that  makes this not work....
+    β = α
+    Cnαβ = 2^(α+β+1)  
+    Cnαβ *= gamma(n + α +  1) * gamma(n + β + 1) / gamma(n + α + β + 1)
+    Cnαβ /= gamma(1 + n)
+    val = Cnαβ / (1-x^2) / dπx^2
+    val
+end
+```
 
 """
 glaser_liu_rokhlin_gauss_nodes()
@@ -122,6 +157,7 @@ function prufer(Π::Type{P},n) where {P <: AbstractCCOP}
         a, b = first(dom)+eps(), last(dom)-eps()
         x = clamp(x, a, b)
         p,q,r,dp,dq,dr = pqr(Π, n, x)
+        @show p,q,r,dp, dq, dr
         -inv(sqrt(r/p) + (dr*p - dp*r + 2r*q)/(2r*p) * sin(2θ)/2)
     end
 end
diff --git a/src/Orthogonal/orthogonal.jl b/src/Orthogonal/orthogonal.jl
index 662da8d..9b45aab 100644
--- a/src/Orthogonal/orthogonal.jl
+++ b/src/Orthogonal/orthogonal.jl
@@ -1,5 +1,7 @@
 ## Abstract  types for  orthogonal  polynomials
 
+export Basis
+
 ## Has An(P), Bn(P), Cn(P)
 abstract type AbstractOrthogonalPolynomial{T} <: AbstractSpecialPolynomial{T} end
 abstract type AbstractContinuousOrthogonalPolynomial{T} <: AbstractOrthogonalPolynomial{T} end
@@ -33,45 +35,31 @@ function Polynomials.fromroots(P::Type{<:AbstractOrthogonalPolynomial}, roots::A
 end
 
 
+Base.convert(P::Type{<:AbstractOrthogonalPolynomial}, c::Number) = c * one(P)
+Base.one(P::Type{<:AbstractOrthogonalPolynomial}, var::Polynomials.SymbolLike=:x) =   basis(P,0, var) / k0(P)
+
 Polynomials.variable(P::Type{<:AbstractOrthogonalPolynomial},  var::Polynomials.SymbolLike=:x) =
-    (basis(P,1,var) - Bn(P,0)) / An(P,0)
+    (basis(P,1,var) / k0(P) - Bn(P,0)) / An(P,0) 
 
 ## Evaluation
 
 # from type, cs, x
 function clenshaw_eval(P::Type{<:AbstractOrthogonalPolynomial{T}}, cs, x::S) where {T,S}
-
     N = length(cs)
-    N == 0 && return zero(T)*zero(S)
-    N == 1 && return cs[1] * one(S)
+    p₀ = k0(P)
+    R = promote_type(promote_type(T,S), typeof(An(P,0)))
+    N == 0 && return zero(R)
+    N == 1 && return (cs[1] * p₀) * one(R) 
 
-    Δ0 = cs[end - 1]
-    Δ1 = cs[end]
+    Δ0::R = cs[end - 1]
+    Δ1::R = cs[end]
     @inbounds for i in N-1:-1:2
         Δ0, Δ1 = cs[i - 1] - Δ1 * Cn(P, i-1), Δ0 + Δ1 * muladd(x, An(P,i-1),Bn(P,i-1))
     end
-
-    return Δ0 + Δ1 * muladd(x, An(P,0),  Bn(P,0))
+    p₁ =  muladd(x, An(P,0),  Bn(P,0)) * p₀
+    return Δ0 * p₀  + Δ1 * p₁
 end
 
-# # from instance, x
-# function clenshaw_eval(p::P, x::S) where {P <: AbstractOrthogonalPolynomial, S}
-
-#     T, cs = eltype(p), coeffs(p)
-#     N = length(cs)
-#     N == 0 && return zero(T)*zero(S)
-#     N == 1 && return cs[1] * one(S)
-
-#     Δ0 = cs[end - 1]
-#     Δ1 = cs[end]
-#     @inbounds for i in N-1:-1:2
-#         Δ0, Δ1 = cs[i - 1] - Δ1 * Cn(P, i-1), Δ0 + Δ1 * muladd(x, An(P,i-1),Bn(P,i-1))
-#     end
-
-#     return Δ0 + Δ1 * muladd(x, An(P,0),  Bn(P,0))
-# end
-
-
 
 ## Connection/Linearization
 
@@ -123,6 +111,7 @@ end
 
 # is P a monic polynomial system?
 ismonic(::Type{P}) where {P <: AbstractOrthogonalPolynomial} = false
+isorthonormal(::Type{P}) where {P <: AbstractOrthogonalPolynomial} = false
 
 # cf. https://en.wikipedia.org/wiki/Orthogonal_polynomials#Recurrence_relation
 # Orthogonal polynomials have a three-point recursion formula
@@ -195,7 +184,29 @@ function monic(p::P) where {P <: AbstractOrthogonalPolynomial}
     p / (p[end]*leading_term(P, n))
 end
 
+##
+## Work with compact  basis
+##
+"""
+     Basis(P,n)
+     Basis{P}(n)
 
+The  command `basis(P,n, [var])` realizes the polynomial. `Basis(P,n)` does  not.  This  can  be  useful for  evaluation.
+"""
+struct Basis{Π}
+    n::Int
+end
+Basis(P,n) =  Basis{P}(n)
+Base.show(io::IO,  mimetype::MIME"text/plain", b::Basis{P})  where  {P} = print(io, "$(P.name)($(b.n))") 
+
+function  innerproduct(::P,  f::Basis{P}, g::Basis{P}) where {P}
+    n,m  = f.n, g.n
+    if n == m
+        return norm2(P, n)
+    else
+        return zero(P)
+    end
+end
 
 
 
@@ -244,13 +255,12 @@ jacobi_matrix(p::P, n) where {P <: AbstractOrthogonalPolynomial} = jacobi_matrix
 
 Returns a tuple of nodes and weights for Gauss quadrature for the given orthogonal type.
 
+When available, the values are computed through  the `FastGaussQuadratures` package.
+
 For some types, a method from  A. Glaser, X. Liu, and V. Rokhlin. "A fast algorithm for the calculation of the roots of special functions." SIAM J. Sci. Comput., 29 (2007), 1420-1438. is used. 
 
 For others the Jacobi matrix, J_n, for which the Golub-Welsch] algorithm The nodes  are computed from the eigenvalues of J_n, the weights a scaling of the first component of the normalized eigen vectors (β_0 * [v[1] for v in vs])
 
-!!! note
-    See the [FastGaussQuadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) package for faster, vastly more engineered implementations.
-
 """
 function gauss_nodes_weights(p::Type{P}, n) where {P <: AbstractOrthogonalPolynomial}
     J = jacobi_matrix(P, n)
@@ -261,6 +271,9 @@ function gauss_nodes_weights(p::Type{P}, n) where {P <: AbstractOrthogonalPolyno
     eig.values,  wts
 end
 
+gauss_nodes_weights(B::Basis{P})  where  {P} = gauss_nodes_weights(B.P, B.n)
+
+    
 
 ##
 ## --------------------------------------------------
diff --git a/src/SpecialPolynomials.jl b/src/SpecialPolynomials.jl
index 9158b08..e01de1d 100644
--- a/src/SpecialPolynomials.jl
+++ b/src/SpecialPolynomials.jl
@@ -11,6 +11,9 @@ export basis
 using QuadGK
 using Memoize
 
+using HypergeometricFunctions
+using FastTransforms
+using FastGaussQuadrature
 
 
 include("utils.jl")
@@ -44,7 +47,6 @@ include("Orthogonal/connection.jl")
 include("Orthogonal/glaser-liu-rokhlin.jl")
 
 
-
 include("Interpolating/interpolating.jl")
 include("Interpolating/Lagrange.jl")
 include("Interpolating/Newton.jl")
diff --git a/src/utils.jl b/src/utils.jl
index ad6b121..96644f4 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -37,6 +37,7 @@ end
 # default is (x)ᵢ
 Pochhammer(z,n::Int) = Pochhammer(Val(:rising),  z, n)
 
+
 # (x)_n/n!
 Pochhammer_factorial(z, n) = Pochhammer_factorial(Val(:rising), z, n)
 
@@ -117,7 +118,17 @@ end
 ## tuples
 @inline plus_1(as) = map(x->x+1, as)
         
+
+# use HypergeometricFunction.mFn if possible
+function pFq(as::Tuple{A,B}, bs::Tuple{C}, z::AbstractFloat) where {A,B,C}
+    HypergeometricFunctions._₂F₁(as[1],as[2],bs[1],z)
+end
+function pFq(as::Tuple{A,B,C}, bs::Tuple{D,E}, z::AbstractFloat) where {A,B,C,D,E}
+    HypergeometricFunctions._₃F₂(as[1],as[2],as[3],bs[1],bs[2],z)
+end
     
+
+
 ## Try to speed up  quadgk by not specializing on F
 mutable struct Wrapper
     F
diff --git a/test/Orthogonal.jl b/test/Orthogonal.jl
index 5c4f4b1..95416a6 100644
--- a/test/Orthogonal.jl
+++ b/test/Orthogonal.jl
@@ -3,17 +3,22 @@
 T = Float64
 Ps = (
     Chebyshev,
+    OrthonormalChebyshev,
     ChebyshevU,
     Laguerre{0},
     Laguerre{1/2},
+    OrthonormalLaguerre{1/2},
     Hermite,
     ChebyshevHermite,
     Jacobi{1/2, 1/2},
     Jacobi{1/4, 3/4},
     Jacobi{-1/2, 1/2},
-    Jacobi{1/2, -1/2},      
+    Jacobi{1/2, -1/2},
+    OrthonormalJacobi{1/2,1/2},
     Legendre,
+    OrthonormalLegendre,
     Gegenbauer{1/2},
+    OrthonormalGegenbauer{1/2},
     Bessel{3/2}, # Bessel{1} is an issue
     Bessel{1/2}
 )
@@ -205,10 +210,30 @@ end
             end
         end
     end
+
+    # Connections via FastTransforms
+    x = variable()
+    for P in Ps
+        for Q in Ps
+            p = P([1.0, 2.0, 3.0])
+            @test convert(Q, p)(x) ≈ p(x)
+        end
+    end
+            
     
 end
 
 
+@testset  "Evaluation" begin
+    for  P in  Ps
+        (SP.ismonic(P) || SP.isorthonormal(P)) && continue
+        n,x = 4, 1.23
+        p = basis(P, n)
+        # compares  clenshaw  to  hypergeometric
+        @test p(x) ≈ Basis(P,n)(x)
+    end
+end
+
 @testset "Arithmetic" begin
 
     x = variable(Polynomial{Float64})
@@ -245,9 +270,9 @@ end
         xs = [1,2,3,4]
         p = P(xs)
         @test -p == P(-xs)
-        ys = copy(xs)
-        ys[1] += 1
-        @test p + 1 == P(ys)
+        ys = copy(xs ./ 1) # might need float
+        ys[1] += one(P)[0]
+        @test p + 1 ≈ p + P(1)
         @test 2p == P(2xs)
     end
 
@@ -318,8 +343,10 @@ end
             ps = rand(1:10, 5)
             p = P(ps)
             @test  (derivative ∘ derivative)(p) ≈ derivative(p, 2)
-            
-            @test (derivative ∘ integrate)(p) ≈ p
+
+            q = (derivative ∘ integrate)(p) - p
+            q = chop(q, atol=1e-8)
+            @test degree(q) <= 0
 
             q = (integrate ∘ derivative)(p) - p
             q = chop(q, atol=1e-8)