Add associated Legendre Polynomials P_l,m (#7)

* Add associated Legendre Polynomials P_l,m

* Correct argument names in documentation

* Correct example output in doc strings

* Correct output in documentation

* New version 0.3.4

Author: eschnett

Project.toml

@@ -1,6 +1,6 @@
 name = "LegendrePolynomials"
 uuid = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
-version = "0.3.3"
+version = "0.3.4"
 
 [deps]
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"

README.md

@@ -26,6 +26,13 @@ julia> Pl(0.5, 20)
 -0.04835838106737356
 ```
 
+To compute the associated Legendre polynomial of degree `l,m` at the argument `x`, use `Plm(x, l, m)`:
+
+```julia
+julia> Plm(0.5, 10, 5)
+30086.169706116176
+```
+
 To compute the n-th derivative of the Legendre polynomial of degree `l` at the argument `x`, use `dnPl(x, l, n)`:
 
 ```julia
@@ -46,6 +53,19 @@ julia> collectPl(0.5, lmax = 5)
   0.08984375
 ```
 
+To compute all the associated Legendre polynomials for `0 <= l <= lmax`, use `collectPlm(x; m, lmax)`
+
+```julia
+julia> collectPlm(0.5, lmax = 5, m = 3)
+6-element OffsetArray(::Vector{Float64}, 0:5) with eltype Float64 with indices 0:5:
+   0.0
+   0.0
+   0.0
+  -9.742785792574933
+ -34.099750274012266
+ -42.62468784251533
+```
+
 To compute all the n-th derivatives for `0 <= l <= lmax`, use `collectdnPl(x; n, lmax)`
 
 ```julia
@@ -63,4 +83,4 @@ Check the documentation for other usage.
 
 # License
 
-This project is licensed under the MIT License - see the [LICENSE](https://github.com/jishnub/LegendrePolynomials.jl/blob/master/LICENSE) file for details.
+This project is licensed under the MIT License - see the [LICENSE](https://github.com/jishnub/LegendrePolynomials.jl/blob/master/LICENSE) file for details.

docs/src/index.md

@@ -7,7 +7,7 @@ end
 
 # Introduction
 
-Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials) and their derivatives using a 3-term recursion relation (Bonnet’s recursion formula).
+Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials), [Associated Legendre polynomials](https://en.wikipedia.org/wiki/Associated_Legendre_polynomials), and their derivatives using a 3-term recursion relation (Bonnet’s recursion formula).
 
 ```math
 P_\ell(x) = \left((2\ell-1) x P_{\ell-1}(x) - (\ell-1)P_{\ell - 2}(x)\right)/\ell
@@ -19,10 +19,12 @@ Currently this package evaluates the standard polynomials that satisfy ``P_\ell(
 \int_{-1}^1 P_m(x) P_n(x) dx = \frac{2}{2n+1} \delta_{mn}.
 ```
 
-There are four main functions:
+There are six main functions: 
 
 * [`Pl(x,l)`](@ref Pl): this evaluates the Legendre polynomial for a given degree `l` at the argument `x`. The argument needs to satisfy `-1 <= x <= 1`.
 * [`collectPl(x; lmax)`](@ref collectPl): this evaluates all the polynomials for `l` lying in `0:lmax` at the argument `x`. As before the argument needs to lie in the domain of validity. Functionally this is equivalent to `Pl.(x, 0:lmax)`, except `collectPl` evaluates the result in one pass, and is therefore faster. There is also the in-place version [`collectPl!`](@ref) that uses a pre-allocated array.
+* [`Plm(x, l, m)`](@ref Plm): this evaluates the associated Legendre polynomial ``P_\ell,m(x)`` at the argument ``x``. The argument needs to satisfy `-1 <= x <= 1`.
+* [`collectPlm(x; m, lmax)`](@ref collectPlm): this evaluates the associated Legendre polynomials with coefficient `m` for `l = 0:lmax`. There is also an in-place version [`collectPlm!`](@ref) that uses a pre-allocated array.
 * [`dnPl(x, l, n)`](@ref dnPl): this evaluates the ``n``-th derivative of the Legendre polynomial ``P_\ell(x)`` at the argument ``x``. The argument needs to satisfy `-1 <= x <= 1`.
 * [`collectdnPl(x; n, lmax)`](@ref collectdnPl): this evaluates the ``n``-th derivative of all the Legendre polynomials for `l = 0:lmax`. There is also an in-place version [`collectdnPl!`](@ref) that uses a pre-allocated array.
 
@@ -35,6 +37,13 @@ julia> Pl(0.5, 3)
 -0.4375
 ```
 
+Evaluate the associated Legendre Polynomial one `l,m` pair as `Plm(x, l, m)`:
+
+```jldoctest
+julia> Plm(0.5, 3, 2)
+5.625
+```
+
 Evaluate the `n`th derivative for one `l` as `dnPl(x, l, n)`:
 
 ```jldoctest
@@ -53,6 +62,19 @@ julia> collectPl(0.5, lmax = 3)
  -0.4375
 ```
 
+Evaluate all the associated Legendre Polynomials for coefficient `m` as `collectPlm(x; lmax, m)`:
+
+```jldoctest
+julia> collectPlm(0.5, lmax = 5, m = 3)
+6-element OffsetArray(::Vector{Float64}, 0:5) with eltype Float64 with indices 0:5:
+   0.0
+   0.0
+   0.0
+  -9.742785792574933
+ -34.099750274012266
+ -42.62468784251533
+```
+
 Evaluate all the `n`th derivatives as `collectdnPl(x; lmax, n)`:
 
 ```jldoctest

src/LegendrePolynomials.jl

@@ -5,6 +5,9 @@ using OffsetArrays
 export Pl
 export collectPl
 export collectPl!
+export Plm
+export collectPlm
+export collectPlm!
 export dnPl
 export collectdnPl
 export collectdnPl!
@@ -27,6 +30,18 @@ end
 	convert(T, Pl)
 end
 
+@inline function Plm_recursion(::Type{T}, ℓ, m, P_m_lm1, P_m_lm2, x) where {T}
+	P_m_l = ((2ℓ-1) * x * P_m_lm1 - (ℓ+m-1) * P_m_lm2) / (ℓ-m)
+	convert(T, P_m_l)
+end
+
+# special case
+# l == m, in which case there are no (l-1,m) and (l-2,m) terms
+@inline function Plm_recursion_m(::Type{T}, ℓ, m, P_mm1_lm1, x) where {T}
+	P_m_l = -(2ℓ-1) * sqrt(1 - x^2) * P_mm1_lm1
+	convert(T, P_m_l)
+end
+
 @inline function dPl_recursion(::Type{T}, ℓ, n, P_n_lm1, P_nm1_lm1, P_n_lm2, x) where {T}
 	P_n_l = ((2ℓ - 1) * (x * P_n_lm1 + n * P_nm1_lm1) - (ℓ - 1) * P_n_lm2)/ℓ
 	convert(T, P_n_l)
@@ -53,7 +68,7 @@ Return an iterator that generates the values of the Legendre polynomials ``P_\\e
 If `lmax` is specified then only the values of ``P_\\ell(x)`` from `0` to `lmax` are returned.
 
 # Examples
-```jldoctest
+``jldoctest
 julia> import LegendrePolynomials: LegendrePolynomialIterator
 
 julia> iter = LegendrePolynomialIterator(0.5, 4);
@@ -69,15 +84,15 @@ julia> collect(iter)
 julia> iter = LegendrePolynomialIterator(0.5);
 
 julia> collect(Iterators.take(iter, 5)) # evaluete 5 elements (l = 0:4)
-5-element Vector{Float64}:
+5-element Array{Float64,1}:
   1.0
   0.5
  -0.125
  -0.4375
  -0.2890625
 
 julia> collect(Iterators.take(Iterators.drop(iter, 100), 5)) # evaluate Pl for l = 100:104
-5-element Vector{Float64}:
+5-element Array{Float64,1}:
  -0.0605180259618612
   0.02196749072249231
   0.08178451892628381
@@ -159,6 +174,92 @@ function doublefactorial(T, n)
 	convert(T, p)
 end
 
+function _checkvalues_m(x, l, m)
+	assertnonnegative(l)
+	checkdomain(x)
+	m >= 0 || throw(ArgumentError("coefficient m must be >= 0"))
+end
+
+Base.@propagate_inbounds function _unsafePlm!(cache, x, l, m)
+	# unsafe, assumes 1-based indexing
+	checklength(cache, l - m + 1)
+	# handle m == 0
+	collectPl!(cache, x, lmax = l - m)
+
+	# handle m > 0
+	for mi in 1:m
+		# We denote the terms as P_mi_li
+	
+		# li == mi
+		P_mim1_mim1 = cache[1]
+		P_mi_mi = Plm_recursion_m(eltype(cache), mi, mi, P_mim1_mim1, x)
+		cache[1] = P_mi_mi
+
+		if l > m
+			# li == mi + 1
+			# P_mim1_mi = cache[2]
+			# P_mi_mip1 = Plm_recursion(eltype(cache), mi + 1, mi, P_mi_mi, P_mim1_mi, nothing, x)
+			# cache[2] = P_mi_mip1
+			P_mi_lim2 = zero(x)
+			P_mi_lim1 = cache[1]
+			P_mi_li = Plm_recursion(eltype(cache), mi + 1, mi, P_mi_lim1, P_mi_lim2, x)
+			cache[2] = P_mi_li
+
+			for li in mi+2:min(l, l - m + mi)
+				P_mi_lim2 = cache[li - mi - 1]
+				P_mi_lim1 = cache[li - mi]
+				P_mi_li = Plm_recursion(eltype(cache), li, mi, P_mi_lim1, P_mi_lim2, x)
+				cache[li - mi + 1] = P_mi_li
+			end
+		end
+	end
+	nothing
+end
+
+"""
+	Plm(x, l::Integer, m::Integer, [cache::AbstractVector])
+
+Compute the associatedLegendre polynomial ``P_\\ell,m(x)``.
+Optionally a pre-allocated vector `cache` may be provided, which must have a minimum length of `l - m + 1` 
+and may be overwritten during the computation.
+
+The coefficient `m` must be non-negative. For `m == 0` this function just returns 
+Legendre polynomials.
+
+# Examples
+
+```jldoctest
+julia> Plm(0.5, 3, 2)
+5.625
+
+julia> Plm(0.5, 4, 0) == Pl(0.5, 4)
+true
+```
+"""
+Base.@propagate_inbounds function Plm(x, l::Integer, m::Integer, 
+	A = begin
+		_checkvalues_m(x, l, m)
+		# do not allocate A if the value is trivially zero
+		if l < m
+			return zero(polytype(x))
+		end
+		zeros(polytype(x), l - m + 1)
+	end
+	)
+
+	_checkvalues_m(x, l, m)
+	# check if the value is trivially zero in case A is provided in the function call
+	if l < m
+		return zero(eltype(A))
+	end
+
+	cache = OffsetArrays.no_offset_view(A)
+	# function barrier, as no_offset_view may be type-unstable
+	_unsafePlm!(cache, x, l, m)
+
+	return cache[l - m + 1]
+end
+
 function _checkvalues(x, l, n)
 	assertnonnegative(l)
 	checkdomain(x)
@@ -304,7 +405,73 @@ julia> collectPl(0.5, lmax = 4)
 collectPl(x; lmax::Integer) = collect(LegendrePolynomialIterator(x, lmax))
 
 """
-	collectdnPl(x; n::Integer, lmax::Integer)
+	collectPlm(x; lmax::Integer, m::Integer)
+
+Compute the associated Legendre Polynomial ``P_\\ell,m(x)`` for the argument `x` and all degrees `l = 0:lmax`. 
+
+The coefficient `m` must be greater than or equal to zero.
+
+Returns `v` with indices `0:lmax`, where `v[l] == Plm(x, l, m)`.
+
+# Examples
+
+```jldoctest
+julia> collectPlm(0.5, lmax = 3, m = 2)
+4-element OffsetArray(::Vector{Float64}, 0:3) with eltype Float64 with indices 0:3:
+ 0.0
+ 0.0
+ 2.25
+ 5.625
+```
+"""
+function collectPlm(x; lmax::Integer, m::Integer)
+	assertnonnegative(lmax)
+	m >= 0 || throw(ArgumentError("coefficient m must be >= 0"))
+	v = zeros(polytype(x), 0:lmax)
+	if lmax >= m
+		collectPlm!(parent(v), x; lmax = lmax, m = m)
+	end
+	v
+end
+
+"""
+	collectPlm!(v::AbstractVector, x; lmax::Integer, m::Integer)
+
+Compute the associated Legendre Polynomial ``P_\\ell,m(x)`` for the argument `x` and all degrees `l = 0:lmax`, 
+and store the result in `v`. 
+
+The coefficient `m` must be greater than or equal to zero.
+
+At output, `v[l + firstindex(v)] == Plm(x, l, m)` for `l = 0:lmax`.
+
+# Examples
+
+```jldoctest
+julia> v = zeros(4);
+
+julia> collectPlm!(v, 0.5, lmax = 3, m = 2)
+4-element Vector{Float64}:
+ 0.0
+ 0.0
+ 2.25
+ 5.625
+```
+"""
+function collectPlm!(v, x; lmax::Integer, m::Integer)
+	assertnonnegative(lmax)
+	m >= 0 || throw(ArgumentError("coefficient m must be >= 0"))
+	checklength(v, lmax + 1)
+	
+	# trivially zero for l < m
+	fill!((@view v[(0:m-1) .+ firstindex(v)]), zero(eltype(v)))
+	# populate the other elements
+	@inbounds Plm(x, lmax, m, @view v[(m:lmax) .+ firstindex(v)])
+	
+	v
+end
+
+"""
+	collectdnPl(x; lmax::Integer, n::Integer)
 
 Compute the ``n``-th derivative of a Legendre Polynomial ``P_\\ell(x)`` for the argument `x` and all degrees `l = 0:lmax`.
 
@@ -345,7 +512,7 @@ At output, `v[l + firstindex(v)] == dnPl(x, l, n)` for `l = 0:lmax`.
 
 # Examples
 
-```jldoctest
+``jldoctest
 julia> v = zeros(4);
 
 julia> collectdnPl!(v, 0.5, lmax = 3, n = 2)