Add tutorial to docs (#14)

* add tutorial to docs

* expand input types

* test with Symbolics

* removed unused polytype method

* fix negative m norm

* update tutorial

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "LegendrePolynomials"
 uuid = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
-version = "0.4.1"
+version = "0.4.2"
 
 [deps]
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
@@ -10,15 +10,19 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Aqua = "0.5"
 HyperDualNumbers = "4"
 OffsetArrays = "1"
+PerformanceTestTools = "0.1"
 QuadGK = "2"
 SpecialFunctions = "1, 2"
+Symbolics = "4"
 julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
+PerformanceTestTools = "dc46b164-d16f-48ec-a853-60448fc869fe"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "HyperDualNumbers", "Test", "QuadGK"]
+test = ["Aqua", "HyperDualNumbers", "PerformanceTestTools", "Symbolics", "Test", "QuadGK"]

docs/Project.toml

@@ -2,7 +2,9 @@
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
 LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 Documenter = "0.27"
 HyperDualNumbers = "4"
+Symbolics = "4"

docs/make.jl

@@ -15,6 +15,7 @@ makedocs(;
     ),
     pages=[
         "Introduction" => "index.md",
+        "Tutorial" => "tutorial.md",
         "Derivatives" => "derivatives.md",
     ],
 )

docs/src/index.md

@@ -48,10 +48,10 @@ P_{\ell}^{m}\left(x\right)=\sqrt{\frac{2\left(\ell+m\right)!}{\left(2\ell+1\righ
 
 There are six main functions:
 
-* [`Pl(x,l; [norm])`](@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, [norm])`](@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; [norm])`](@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, [norm])`](@ref collectPlm): this evaluates the associated Legendre polynomials with coefficient `m` for `l = abs(m):lmax`. There is also an in-place version [`collectPlm!`](@ref) that uses a pre-allocated array.
+* [`Pl(x,l; [norm = Val(:standard)])`](@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, [lmin = 0], [norm = Val(:standard)])`](@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, lmin: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; [norm = Val(:standard)])`](@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, [lmin = abs(m)], [norm = Val(:standard)])`](@ref collectPlm): this evaluates the associated Legendre polynomials with coefficient `m` for `l = lmin: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.
 
@@ -77,7 +77,7 @@ julia> p ≈ 45/8 # analytical value
 true
 ```
 
-Evaluate the `n`th derivative for one `l` as `dnPl(x, l, n)`:
+Evaluate the `n`th derivative of `Pl` for one `l` as `dnPl(x, l, n)`:
 
 ```jldoctest
 julia> dnPl(0.5, 3, 2)
@@ -107,7 +107,7 @@ julia> collectPlm(0.5, lmax = 5, m = 3)
  -42.62468784251535
 ```
 
-Evaluate all the `n`th derivatives as `collectdnPl(x; lmax, n)`:
+Evaluate all the `n`th derivatives of `Pl` as `collectdnPl(x; lmax, n)`:
 
 ```jldoctest
 julia> collectdnPl(0.5, lmax = 5, n = 3)
@@ -120,43 +120,6 @@ julia> collectdnPl(0.5, lmax = 5, n = 3)
  65.625
 ```
 
-# Increase precision
-
-The precision of the result may be changed by using arbitrary-precision types such as `BigFloat`. For example, using `Float64` arguments we obtain
-
-```jldoctest
-julia> Pl(1/3, 5)
-0.33333333333333337
-```
-
-whereas using `BigFloat`, we obtain
-
-```jldoctest
-julia> Pl(big(1)/3, 5)
-0.3333333333333333333333333333333333333333333333333333333333333333333333333333305
-```
-
-The precision of the latter may be altered using `setprecision`, as
-
-```jldoctest
-julia> setprecision(300) do
-       Pl(big(1)/3, 5)
-       end
-0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333317
-```
-
-This is particularly important to avoid overflow while computing high-order derivatives. For example:
-
-```jldoctest
-julia> dnPl(0.5, 300, 200) # Float64
-NaN
-
-julia> dnPl(big(1)/2, 300, 200) # BigFloat
-1.738632750542319394663553898425873258768856732308227932150592526951212145232716e+499
-```
-
-In general, one would need to use higher precision for both the argument `x` and the degree `l` to obtain accurate results.
-
 # Reference
 
 ```@autodocs

docs/src/tutorial.md

@@ -0,0 +1,137 @@
+```@meta
+DocTestSetup = :(using LegendrePolynomials)
+```
+
+## Computing normalized Legendre polynomials
+
+By default, the Legendre and associated Legendre polynomials are not normalized.
+One may specify the normalization through the keyword argument `norm`.
+The normalization options accepted are
+
+* `norm = Val(:standard)`: standard, unnormalized polynomials. This is the default option.
+* `norm = Val(:normalized)`: fully normalized polynomials with an L2 norm of 1
+
+!!! note
+	Irrespective of the norm specified, the 3-term recursion relations used are stable ones,
+	so polynomials of high degrees may be computed without overflow.
+
+```jldoctest
+julia> l = 3;
+
+julia> Pl(0.5, l)
+-0.4375
+
+julia> Pl(0.5, l, norm = Val(:normalized))
+-0.8184875533567997
+
+julia> Pl(0.5, l, norm = Val(:normalized)) ≈ √((2l+1)/2) * Pl(0.5, l)
+true
+
+julia> l = m = 3000;
+
+julia> Plm(0.5, l, m, norm = Val(:normalized))
+2.172276347346834e-187
+```
+
+## Condon-Shortley phase
+
+By default, the associated Legendre polynomials are computed by including the Condon-Shortley phase ``(-1)^m``.
+This may be disabled by passing the flag `csphase = false` to the various `Plm` functions.
+
+!!! note
+    The `Pl` functions do not accept the keyword argument `csphase`.
+
+```jldoctest
+julia> Plm(0.5, 3, 3)
+-9.742785792574933
+
+julia> Plm(0.5, 3, 3, csphase = false)
+9.742785792574933
+
+julia> all(Plm(0.5, 3, m, csphase = false) ≈ (-1)^m * Plm(0.5, 3, m) for m in -3:3)
+true
+```
+
+## Polynomials at multiple points
+
+The `collectPl(m)` functions return a vector of polynomials evaluated for one ``\theta``. To evaluate the polynomials for multiple ``\theta`` in one go, one may run
+
+```jldoctest multipletheta
+julia> θr = range(0, pi, length=6);
+
+julia> collectPl.(cos.(θr), lmax=3)
+6-element Vector{OffsetArrays.OffsetVector{Float64, Vector{Float64}}}:
+ [1.0, 1.0, 1.0, 1.0]
+ [1.0, 0.8090169943749475, 0.4817627457812106, 0.1102457514062632]
+ [1.0, 0.30901699437494745, -0.35676274578121053, -0.38975424859373686]
+ [1.0, -0.30901699437494734, -0.35676274578121064, 0.3897542485937368]
+ [1.0, -0.8090169943749473, 0.48176274578121037, -0.11024575140626285]
+ [1.0, -1.0, 1.0, -1.0]
+```
+
+This returns a vector of vectors, where each element corresponds to one ``\theta``. Often, one wants a `Matrix` where
+each column corresponds to one ``\theta``. We may obtain this as
+
+```jldoctest multipletheta
+julia> mapreduce(hcat, θr) do θ
+           collectPl(cos(θ), lmax=3)
+       end
+4×6 Matrix{Float64}:
+ 1.0  1.0        1.0        1.0        1.0        1.0
+ 1.0  0.809017   0.309017  -0.309017  -0.809017  -1.0
+ 1.0  0.481763  -0.356763  -0.356763   0.481763   1.0
+ 1.0  0.110246  -0.389754   0.389754  -0.110246  -1.0
+```
+
+As of Julia `v1.7.2` and `OffsetArrays` `v1.10.8`, this strips off the axis information, so one would need to wrap the result in an `OffsetArray` again to index the `Matrix` using the degrees ``\ell``.
+
+## Increasing precision
+
+The precision of the result may be changed by using arbitrary-precision types such as `BigFloat`. For example, using `Float64` arguments we obtain
+
+```jldoctest
+julia> Pl(1/3, 5)
+0.33333333333333337
+```
+
+whereas using `BigFloat`, we obtain
+
+```jldoctest
+julia> Pl(big(1)/3, 5)
+0.3333333333333333333333333333333333333333333333333333333333333333333333333333305
+```
+
+This is particularly important to avoid overflow while computing high-order derivatives. For example:
+
+```jldoctest
+julia> dnPl(0.5, 300, 200) # Float64
+NaN
+
+julia> dnPl(big(1)/2, 300, 200) # BigFloat
+1.738632750542319394663553898425873258768856732308227932150592526951212145232716e+499
+```
+
+In general, one would need to use higher precision for both the argument `x` and the degrees `l` to obtain accurate results. For example:
+
+```jldoctest
+julia> Plm(big(1)/2, big(3000), 3000)
+2.05451584347939475644802290993338963448971107938391335496027846832433343889916e+9844
+```
+
+## Symbolic evaluation
+
+It's possible to symbolically evaluate Legendre or associated Legendre polynomials using [`Symbolics.jl`](https://github.com/JuliaSymbolics/Symbolics.jl). An exmaple is:
+
+```jldoctest
+julia> using Symbolics
+
+julia> @variables x;
+
+julia> Pl(x, 3)
+(5//3)*x*((3//2)*(x^2) - (1//2)) - (2//3)*x
+
+julia> myf = eval(build_function(Pl(x, 3), [x]));
+
+julia> myf(0.4) == Pl(0.4, 3)
+true
+```

src/LegendrePolynomials.jl

@@ -28,9 +28,14 @@ function __init__()
     end
 end
 
-function checkdomain(x)
-    abs(x) > 1 && throw(DomainError(x,
+checkdomain(x) = true
+function checkdomain(x::Number)
+    f = abs(x) > 1
+    if f isa Bool
+        f && throw(DomainError(x,
         "Legendre Polynomials are defined for arguments lying in -1 ⩽ x ⩽ 1"))
+    end
+    return nothing
 end
 
 assertnonnegative(l::Integer) = l >= 0 || throw(ArgumentError("l must be >= 0, received " * string(l)))
@@ -119,9 +124,13 @@ function coeff(l, m, T = float(promote_type(typeof(l), typeof(m))))
     √((2T(l) - 1)/((T(l)-m)*(T(l)+m))*(2T(l)+1))
 end
 
-polytype(x) = float(typeof(x))
-polytype(m::Nothing) = Float64
-polytype(x, m) = promote_type(polytype(x), polytype(m))
+polytype(x) = typeof(float.(x))
+polytype(x::Array{<:Number}) = Array{promote_type(eltype(x), Float64), ndims(x)}
+polytype(x::Number) = float(typeof(x))
+polytype(x::Number, m::Number) = promote_type(polytype(x), polytype(m))
+polytype(x, m::Number) = typeof(float.(x) .* m)
+polytype(x::Array{<:Number}, m::Number) = Array{promote_type(eltype(x), Float64, typeof(m)), ndims(x)}
+polytype(x, ::Nothing) = polytype(x)
 
 struct LegendrePolynomialIterator{T, M<:Union{Integer, Nothing}, V}
     x :: V
@@ -145,8 +154,8 @@ Base.IteratorSize(::Type{<:LegendrePolynomialIterator}) = Base.IsInfinite()
 # Iteration for Legendre polynomials
 # starting value, l == 0 for m == 0
 function Base.iterate(iter::LegendrePolynomialIterator{T,Nothing}) where {T}
-    Pl = one(T)
-    Plm1 = zero(T)
+    Pl = convert(T, one(iter.x))
+    Plm1 = convert(T, zero(iter.x))
     nextstate = (Pl, Plm1, nothing)
     return Pl, (1, nextstate)
 end
@@ -166,13 +175,13 @@ function Base.iterate(iter::LegendrePolynomialIterator{T,<:Integer}) where {T}
     m = iter.m
     x = iter.x
     if m == 0
-        Pl = one(T)
+        Pl = convert(T, one(x))
     else
         f = pll_prefactor(m, csphase = iter.csphase)
-        t = f * (√(1-x^2))^m
+        t = f * (√(one(x) - x^2))^m
         Pl = convert(T, t)
     end
-    Plm1 = zero(T)
+    Plm1 = convert(T, zero(x))
     α_ℓ_m = Inf
     return Pl, (m+1, (Pl, Plm1, α_ℓ_m))
 end
@@ -206,7 +215,7 @@ function maybenormalize(P, l, m, norm::Val{:standard}, csphase = true)
     if m == 0
         return maybenormalize(P, l, norm)
     else
-        return (m < 0 && csphase ? neg1pow(m) : 1) * plm_norm(l, m) * P
+        return (m < 0 ? neg1pow(m) : 1) * plm_norm(l, m) * P
     end
 end
 maybenormalize(P, l, m, norm, args...) = throw(ArgumentError("norm = $norm undefined, valid norms are Val(:standard) and Val(:normalized)"))

test/Project.toml

@@ -1,12 +1,17 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+PerformanceTestTools = "dc46b164-d16f-48ec-a853-60448fc869fe"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Aqua = "0.5"
 HyperDualNumbers = "4"
+PerformanceTestTools = "0.1"
 QuadGK = "2"
+Symbolics = "4"