Update Interface for modified Chebyshev moment-based quadrature

For the two strings `clenshawcurtis` and `fejer`, there are now three methods each:

`plan_*`, `*nodes`, `*weights`, where the Fejer methods are followed by a `1` or `2` to indicate the type.

This is a better interface because `*nodes` is more extensible: different Chebyshev moments may be added for new quadrature rules. In particular, this commit also adds `log` moments.

Author: MikaelSlevinsky

docs/src/index.md

@@ -134,6 +134,8 @@ sphevaluate
 
 ## Internal Methods
 
+### Miscellaneous Special Functions
+
 ```@docs
 FastTransforms.half
 ```
@@ -162,38 +164,58 @@ FastTransforms.pochhammer
 FastTransforms.stirlingseries
 ```
 
+### Modified Chebyshev Moment-Based Quadrature
+
 ```@docs
-FastTransforms.clenshawcurtis
+FastTransforms.clenshawcurtisnodes
 ```
 
 ```@docs
 FastTransforms.clenshawcurtisweights
 ```
 
 ```@docs
-FastTransforms.fejer1
+FastTransforms.fejernodes1
 ```
 
 ```@docs
-FastTransforms.fejer2
+FastTransforms.fejerweights1
 ```
 
 ```@docs
-FastTransforms.fejerweights1
+FastTransforms.fejernodes2
 ```
 
 ```@docs
 FastTransforms.fejerweights2
 ```
 
+```@docs
+FastTransforms.chebyshevmoments1
+```
+
 ```@docs
 FastTransforms.chebyshevjacobimoments1
 ```
 
+```@docs
+FastTransforms.chebyshevlogmoments1
+```
+
+```@docs
+FastTransforms.chebyshevmoments2
+```
+
 ```@docs
 FastTransforms.chebyshevjacobimoments2
 ```
 
+```@docs
+FastTransforms.chebyshevlogmoments2
+```
+
+### Jacobi Polynomial Increment and Decrement Operators
+
 ```@docs
 FastTransforms.incrementα!
 ```

src/ChebyshevJacobiPlan.jl

@@ -175,7 +175,7 @@ function BackwardChebyshevJacobiPlan{T}(c_cheb::AbstractVector{T},α::T,β::T,M:
 
     # Clenshaw-Curtis nodes and weights
     θ = N > 0 ? T[k/2N for k=zero(T):2N] : T[0]
-    w = N > 0 ? clenshawcurtisweights(2N+1,α,β,p₁) : T[0]
+    w = N > 0 ? clenshawcurtisweights!(chebyshevjacobimoments1(T, 2N+1, α, β), p₁) : T[0]
 
     # Initialize sines and cosines
     tempsin = sinpi.(θ./2)

src/ChebyshevUltrasphericalPlan.jl

@@ -169,7 +169,7 @@ function BackwardChebyshevUltrasphericalPlan{T}(c_ultra::AbstractVector{T},λ::T
 
     # Clenshaw-Curtis nodes and weights
     θ = N > 0 ? T[k/2N for k=zero(T):2N] : T[0]
-    w = N > 0 ? clenshawcurtisweights(2N+1,λ-half(λ),λ-half(λ),p₁) : T[0]
+    w = N > 0 ? clenshawcurtisweights!(chebyshevjacobimoments1(T, 2N+1, λ-half(λ), λ-half(λ)), p₁) : T[0]
 
     # Initialize sines and cosines
     tempsin = sinpi.(θ)

src/FastTransforms.jl

@@ -51,9 +51,9 @@ export triones, trizeros, trirand, trirandn, trievaluate
 #export ChebyshevJacobiPlan, jac2cheb, cheb2jac
 #export sqrtpi, pochhammer, stirlingseries, stirlingremainder, Aratio, Cratio, Anαβ
 #export Cnmαβ, Cnαβ, Cnmλ, Cnλ, Λ, absf, findmindices!
-#export clenshawcurtis, clenshawcurtis_plan, clenshawcurtisweights
-#export fejer1, fejer_plan1, fejerweights1
-#export fejer2, fejer_plan2, fejerweights2
+#export plan_clenshawcurtis, clenshawcurtisnodes, clenshawcurtisweights
+#export plan_fejer1, fejernodes1, fejerweights1
+#export plan_fejer2, fejernodes2, fejerweights2
 #export RecurrencePlan, forward_recurrence!, backward_recurrence
 
 include("stepthreading.jl")

src/clenshawcurtis.jl

@@ -1,20 +1,20 @@
-clenshawcurtis_plan(μ) = length(μ) > 1 ? FFTW.plan_r2r!(μ, FFTW.REDFT00) : ones(μ)'
+plan_clenshawcurtis(μ) = length(μ) > 1 ? FFTW.plan_r2r!(μ, FFTW.REDFT00) : ones(μ)'
 
 doc"""
-Compute nodes and weights of the Clenshaw—Curtis quadrature rule with a Jacobi weight.
+Compute nodes of the Clenshaw—Curtis quadrature rule.
 """
-clenshawcurtis{T<:AbstractFloat}(N::Int,α::T,β::T) = clenshawcurtis(N,α,β,clenshawcurtis_plan(zeros(T,N)))
-clenshawcurtis{T<:AbstractFloat}(N::Int,α::T,β::T,plan) = T[cospi(k/(N-one(T))) for k=0:N-1],clenshawcurtisweights(N,α,β,plan)
+clenshawcurtisnodes(::Type{T}, N::Int) where T = T[cospi(k/(N-one(T))) for k=0:N-1]
 
 doc"""
-Compute weights of the Clenshaw—Curtis quadrature rule with a Jacobi weight.
+Compute weights of the Clenshaw—Curtis quadrature rule with modified Chebyshev moments of the first kind ``\mu``.
 """
-clenshawcurtisweights{T<:AbstractFloat}(N::Int,α::T,β::T) = clenshawcurtisweights(N,α,β,clenshawcurtis_plan(zeros(T,N)))
-function clenshawcurtisweights{T<:AbstractFloat}(N::Int,α::T,β::T,plan)
-    μ = chebyshevjacobimoments1(N,α,β)
-    scale!(μ,inv(N-one(T)))
+clenshawcurtisweights(μ::Vector) = clenshawcurtisweights!(copy(μ))
+clenshawcurtisweights!(μ::Vector) = clenshawcurtisweights!(μ, plan_clenshawcurtis(μ))
+function clenshawcurtisweights!(μ::Vector{T}, plan) where T
+    N = length(μ)
+    scale!(μ, inv(N-one(T)))
     plan*μ
-    μ[1]/=2;μ[N]/=2
+    μ[1] *= half(T); μ[N] *= half(T)
     return μ
 end
 

src/fejer.jl

@@ -1,40 +1,38 @@
-fejer_plan1(μ) = FFTW.plan_r2r!(μ, FFTW.REDFT01)
-fejer_plan2(μ) = FFTW.plan_r2r!(μ, FFTW.RODFT00)
+plan_fejer1(μ) = FFTW.plan_r2r!(μ, FFTW.REDFT01)
 
 doc"""
-Compute nodes and weights of Fejer's first quadrature rule with a Jacobi weight.
+Compute nodes of Fejer's first quadrature rule.
 """
-fejer1{T<:AbstractFloat}(N::Int,α::T,β::T) = fejer1(N,α,β,fejer_plan1(zeros(T,N)))
+fejernodes1(::Type{T}, N::Int) where T = T[sinpi((N-2k-one(T))/2N) for k=0:N-1]
 
 doc"""
-Compute nodes and weights of Fejer's second quadrature rule with a Jacobi weight.
+Compute weights of Fejer's first quadrature rule with modified Chebyshev moments of the first kind ``\mu``.
 """
-fejer2{T<:AbstractFloat}(N::Int,α::T,β::T) = fejer2(N,α,β,fejer_plan2(zeros(T,N)))
+fejerweights1(μ::Vector) = fejerweights1!(copy(μ))
+fejerweights1!(μ::Vector) = fejerweights1!(μ, plan_fejer1(μ))
+function fejerweights1!(μ::Vector{T}, plan) where T
+    N = length(μ)
+    scale!(μ, inv(T(N)))
+    return plan*μ
+end
 
-fejer1{T<:AbstractFloat}(N::Int,α::T,β::T,plan) = T[sinpi((N-2k-one(T))/2N) for k=0:N-1],fejerweights1(N,α,β,plan)
-fejer2{T<:AbstractFloat}(N::Int,α::T,β::T,plan) = T[cospi((k+one(T))/(N+one(T))) for k=0:N-1],fejerweights2(N,α,β,plan)
 
+plan_fejer2(μ) = FFTW.plan_r2r!(μ, FFTW.RODFT00)
 
 doc"""
-Compute weights of Fejer's first quadrature rule with a Jacobi weight.
+Compute nodes of Fejer's second quadrature rule.
 """
-fejerweights1{T<:AbstractFloat}(N::Int,α::T,β::T) = fejerweights1(N,α,β,fejer_plan1(zeros(T,N)))
+fejernodes2(::Type{T}, N::Int) where T = T[cospi((k+one(T))/(N+one(T))) for k=0:N-1]
 
 doc"""
-Compute weights of Fejer's second quadrature rule with a Jacobi weight.
+Compute weights of Fejer's second quadrature rule with modified Chebyshev moments of the second kind ``\mu``.
 """
-fejerweights2{T<:AbstractFloat}(N::Int,α::T,β::T) = fejerweights2(N,α,β,fejer_plan2(zeros(T,N)))
-
-function fejerweights1{T<:AbstractFloat}(N::Int,α::T,β::T,plan)
-    μ = chebyshevjacobimoments1(N,α,β)
-    scale!(μ,inv(T(N)))
-    return plan*μ
-end
-
-function fejerweights2{T<:AbstractFloat}(N::Int,α::T,β::T,plan)
-    μ = chebyshevjacobimoments2(N,α,β)
+fejerweights2(μ::Vector) = fejerweights2!(copy(μ))
+fejerweights2!(μ::Vector) = fejerweights2!(μ, plan_fejer2(μ))
+function fejerweights2!(μ::Vector{T}, plan) where T
+    N = length(μ)
     Np1 = N+one(T)
-    scale!(μ,inv(Np1))
+    scale!(μ, inv(Np1))
     plan*μ
     @inbounds for i=1:N μ[i] = sinpi(i/Np1)*μ[i] end
     return μ

src/specialfunctions.jl

@@ -404,15 +404,30 @@ function init_c₁c₂!(c₁::Vector,c₂::Vector,u::Vector,v::Vector,c::Vector,
 end
 
 
+doc"""
+Modified Chebyshev moments of the first kind:
+
+```math
+    \int_{-1}^{+1} T_n(x) {\rm\,d}x.
+```
+"""
+function chebyshevmoments1(::Type{T}, N::Int) where T
+    μ = zeros(T, N)
+    for i = 0:2:N-1
+        @inbounds μ[i+1] = two(T)/T(1-i^2)
+    end
+    μ
+end
+
 doc"""
 Modified Chebyshev moments of the first kind with respect to the Jacobi weight:
 
 ```math
     \int_{-1}^{+1} T_n(x) (1-x)^\alpha(1+x)^\beta{\rm\,d}x.
 ```
 """
-function chebyshevjacobimoments1{T<:AbstractFloat}(N::Int,α::T,β::T)
-    μ = zeros(T,N)
+function chebyshevjacobimoments1(::Type{T}, N::Int, α, β) where T
+    μ = zeros(T, N)
     N > 0 && (μ[1] = 2 .^ (α+β+1)*beta(α+1,β+1))
     if N > 1
         μ[2] = μ[1]*(β-α)/(α+β+2)
@@ -423,15 +438,50 @@ function chebyshevjacobimoments1{T<:AbstractFloat}(N::Int,α::T,β::T)
     μ
 end
 
+doc"""
+Modified Chebyshev moments of the first kind with respect to the logarithmic weight:
+
+```math
+    \int_{-1}^{+1} T_n(x) \log\left(\frac{1-x}{2}\right){\rm\,d}x.
+```
+"""
+function chebyshevlogmoments1(::Type{T}, N::Int) where T
+    μ = zeros(T, N)
+    N > 0 && (μ[1] = -two(T))
+    if N > 1
+        μ[2] = -one(T)
+        for i=1:N-2
+            cst = isodd(i) ? T(4)/T(i^2-4) : T(4)/T(i^2-1)
+            @inbounds μ[i+2] = ((i-2)*μ[i]+cst)/(i+2)
+        end
+    end
+    μ
+end
+
+doc"""
+Modified Chebyshev moments of the second kind:
+
+```math
+    \int_{-1}^{+1} U_n(x) {\rm\,d}x.
+```
+"""
+function chebyshevmoments2(::Type{T}, N::Int) where T
+    μ = zeros(T, N)
+    for i = 0:2:N-1
+        @inbounds μ[i+1] = two(T)/T(i+1)
+    end
+    μ
+end
+
 doc"""
 Modified Chebyshev moments of the second kind with respect to the Jacobi weight:
 
 ```math
     \int_{-1}^{+1} U_n(x) (1-x)^\alpha(1+x)^\beta{\rm\,d}x.
 ```
 """
-function chebyshevjacobimoments2{T<:AbstractFloat}(N::Int,α::T,β::T)
-    μ = zeros(T,N)
+function chebyshevjacobimoments2(::Type{T}, N::Int, α, β) where T
+    μ = zeros(T, N)
     N > 0 && (μ[1] = 2 .^ (α+β+1)*beta(α+1,β+1))
     if N > 1
         μ[2] = 2μ[1]*(β-α)/(α+β+2)
@@ -442,6 +492,24 @@ function chebyshevjacobimoments2{T<:AbstractFloat}(N::Int,α::T,β::T)
     μ
 end
 
+doc"""
+Modified Chebyshev moments of the second kind with respect to the logarithmic weight:
+
+```math
+    \int_{-1}^{+1} U_n(x) \log\left(\frac{1-x}{2}\right){\rm\,d}x.
+```
+"""
+function chebyshevlogmoments2(::Type{T}, N::Int) where T
+    μ = chebyshevlogmoments1(T, N)
+    if N > 1
+        μ[2] *= two(T)
+        for i=1:N-2
+            @inbounds μ[i+2] = 2μ[i+2] + μ[i]
+        end
+    end
+    μ
+end
+
 doc"""
 Compute Jacobi expansion coefficients in ``P_n^{(\alpha+1,\beta)}(x)`` given Jacobi expansion coefficients in ``P_n^{(\alpha,\beta)}(x)`` in-place.
 """

test/basictests.jl

@@ -1,6 +1,8 @@
 using Compat, FastTransforms, LowRankApprox
 using Compat.Test
-import FastTransforms: Cnλ, Λ, lambertw, Cnαβ, Anαβ, fejer1, fejer2, clenshawcurtis
+import FastTransforms: Cnλ, Λ, lambertw, Cnαβ, Anαβ
+import FastTransforms: clenshawcurtisnodes, clenshawcurtisweights, fejernodes1, fejerweights1, fejernodes2, fejerweights2
+import FastTransforms: chebyshevmoments1, chebyshevmoments2, chebyshevjacobimoments1, chebyshevjacobimoments2, chebyshevlogmoments1, chebyshevlogmoments2
 
 @testset "Special functions" begin
     n = 0:1000_000
@@ -33,19 +35,44 @@ end
     N = 20
     f(x) = exp(x)
 
-    x,w = fejer1(N,0.,0.)
-    @test norm(dot(f.(x),w)-2sinh(1)) ≤ 4eps()
-    x,w = fejer2(N,0.,0.)
-    @test norm(dot(f.(x),w)-2sinh(1)) ≤ 4eps()
-    x,w = clenshawcurtis(N,0.,0.)
+    x = clenshawcurtisnodes(Float64, N)
+    μ = chebyshevmoments1(Float64, N)
+    w = clenshawcurtisweights(μ)
     @test norm(dot(f.(x),w)-2sinh(1)) ≤ 4eps()
 
-    x,w = fejer1(N,0.25,0.35)
+    μ = chebyshevjacobimoments1(Float64, N, 0.25, 0.35)
+    w = clenshawcurtisweights(μ)
     @test norm(dot(f.(x),w)-2.0351088204147243) ≤ 4eps()
-    x,w = fejer2(N,0.25,0.35)
+
+    μ = chebyshevlogmoments1(Float64, N)
+    w = clenshawcurtisweights(μ)
+    @test norm(sum(w./(x-3)) - π^2/12) ≤ 4eps()
+
+    x = fejernodes1(Float64, N)
+    μ = chebyshevmoments1(Float64, N)
+    w = fejerweights1(μ)
+    @test norm(dot(f.(x),w)-2sinh(1)) ≤ 4eps()
+
+    μ = chebyshevjacobimoments1(Float64, N, 0.25, 0.35)
+    w = fejerweights1(μ)
     @test norm(dot(f.(x),w)-2.0351088204147243) ≤ 4eps()
-    x,w = clenshawcurtis(N,0.25,0.35)
+
+    μ = chebyshevlogmoments1(Float64, N)
+    w = fejerweights1(μ)
+    @test norm(sum(w./(x-3)) - π^2/12) ≤ 4eps()
+
+    x = fejernodes2(Float64, N)
+    μ = chebyshevmoments2(Float64, N)
+    w = fejerweights2(μ)
+    @test norm(dot(f.(x),w)-2sinh(1)) ≤ 4eps()
+
+    μ = chebyshevjacobimoments2(Float64, N, 0.25, 0.35)
+    w = fejerweights2(μ)
     @test norm(dot(f.(x),w)-2.0351088204147243) ≤ 4eps()
+
+    μ = chebyshevlogmoments2(Float64, N)
+    w = fejerweights2(μ)
+    @test norm(sum(w./(x-3)) - π^2/12) ≤ 4eps()
 end
 
 @testset "Allocation-free ID matrix-vector products" begin