Accept normalization options (#10)

* normalized legendre polynomials

* separate iterations for Pl and Plm

* bugfix in Pl, fix docs

* update docs

* remove unused function normtype

* remove lmax from iterator

* accept lmin in collectPl* functions

* embarrassing bugfix in Pl norm

* Add logfactorial caches

* Add tests for norm

* update docs

* Optionally set the CS phase

* Add test for non-overflow

* promote lmin/lmax in collectPl*

* Add test for large-degree Pl

* Accept negative orders

* Add norm tests for negative order

* update docs for negative m

* Add precision tests for BigFloat

* Add tests for parity

Author: jishnub

.github/workflows/ci.yml

@@ -23,7 +23,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1'
           - 'nightly'
         os:

Project.toml

@@ -1,20 +1,24 @@
 name = "LegendrePolynomials"
 uuid = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
-version = "0.3.6"
+version = "0.4.0"
 
 [deps]
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 Aqua = "0.5"
 HyperDualNumbers = "4"
-OffsetArrays = "0.11, 1"
-julia = "1"
+OffsetArrays = "1"
+QuadGK = "2"
+SpecialFunctions = "2.1"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "HyperDualNumbers", "Test"]
+test = ["Aqua", "HyperDualNumbers", "Test", "QuadGK"]

README.md

@@ -11,7 +11,7 @@ Iterative computation of Legendre Polynomials
 
 ## Installing
 
-To install the package, run 
+To install the package, run
 
 ```julia
 ] add LegendrePolynomials
@@ -53,17 +53,14 @@ julia> collectPl(0.5, lmax = 5)
   0.08984375
 ```
 
-To compute all the associated Legendre polynomials for `0 <= l <= lmax`, use `collectPlm(x; m, lmax)`
+To compute all the associated Legendre polynomials for `abs(m) <= 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
+3-element OffsetArray(::Vector{Float64}, 3:5) with eltype Float64 with indices 3:5:
   -9.742785792574933
- -34.099750274012266
- -42.62468784251533
+ -34.09975027401223
+ -42.62468784251535
 ```
 
 To compute all the n-th derivatives for `0 <= l <= lmax`, use `collectdnPl(x; n, lmax)`

docs/src/index.md

@@ -13,18 +13,45 @@ Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomial
 P_\ell(x) = \left((2\ell-1) x P_{\ell-1}(x) - (\ell-1)P_{\ell - 2}(x)\right)/\ell
 ```
 
-Currently this package evaluates the standard polynomials that satisfy ``P_\ell(1) = 1`` and ``P_0(x) = 1``. These are normalized as
+By default this package evaluates the standard polynomials that satisfy ``P_\ell(1) = 1`` and ``P_0(x) = 1``. These are normalized as
 
 ```math
 \int_{-1}^1 P_m(x) P_n(x) dx = \frac{2}{2n+1} \delta_{mn}.
 ```
 
-There are six main functions: 
+Optionally, normalized polynomials may be evaluated that have an L2 norm of `1`.
 
-* [`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.
+Analogous to Legendre polynomials, one may evaluate associated Legendre polynomials using a 3-term recursion relation. This is evaluated by iterating over the normalized associated Legendre functions, and multiplying the norm at the final stage. Such an iteration avoids floating-point overflow.
+
+The relation used in evaluating normalized associated Legendre polynomials ``\bar{P}_{\ell}^{m}\left(x\right)`` is
+
+```math
+\bar{P}_{\ell}^{m}\left(x\right)=\alpha_{\ell m}\left(x\bar{P}_{\ell-1}^{m}\left(x\right)-\frac{1}{\alpha_{\ell-1m}}\bar{P}_{\ell-2}^{m}\left(x\right)\right),
+```
+
+where
+```math
+\alpha_{\ell m}=\sqrt{\frac{\left(2\ell+1\right)\left(2\ell-1\right)}{\left(\ell-m\right)\left(\ell+m\right)}},
+```
+
+The starting value for a specific ``m`` is obtained by setting ``\ell = m``, in which case we use the explicit relation
+
+```math
+\bar{P}_{m}^{m}\left(x\right)=\left(-1\right)^{m}\left(2m-1\right)!!\sqrt{\frac{\left(2m+1\right)}{2\left(2m\right)!}}\left(1-x^{2}\right)^{\left(m/2\right)}.
+```
+
+The polynomials, thus defined, include the Condon-Shortley phase ``(-1)^m``, and may be evaluated using the standard normalization as
+
+```math
+P_{\ell}^{m}\left(x\right)=\sqrt{\frac{2\left(\ell+m\right)!}{\left(2\ell+1\right)\left(\ell-m\right)!}}\bar{P}_{\ell}^{m}\left(x\right).
+```
+
+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.
 * [`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.
 
@@ -33,15 +60,21 @@ There are six main functions:
 Evaluate the Legendre polynomial for one `l` at an argument`x` as `Pl(x, l)`:
 
 ```jldoctest
-julia> Pl(0.5, 3)
+julia> p = Pl(0.5, 3)
 -0.4375
+
+julia> p ≈ -7/16 # analytical value
+true
 ```
 
 Evaluate the associated Legendre Polynomial one `l,m` pair as `Plm(x, l, m)`:
 
 ```jldoctest
-julia> Plm(0.5, 3, 2)
-5.625
+julia> p = Plm(0.5, 3, 2)
+5.624999999999997
+
+julia> p ≈ 45/8 # analytical value
+true
 ```
 
 Evaluate the `n`th derivative for one `l` as `dnPl(x, l, n)`:
@@ -51,7 +84,8 @@ julia> dnPl(0.5, 3, 2)
 7.5
 ```
 
-Evaluate all the polynomials for `l` in `0:lmax` as `collectPl(x; lmax)`
+Evaluate the Legendre polynomials for `l` in `lmin:lmax` as `collectPl(x; lmin, lmax)`.
+By default `lmin` is chosen to be `0`, and may be omitted.
 
 ```jldoctest
 julia> collectPl(0.5, lmax = 3)
@@ -62,17 +96,15 @@ julia> collectPl(0.5, lmax = 3)
  -0.4375
 ```
 
-Evaluate all the associated Legendre Polynomials for coefficient `m` as `collectPlm(x; lmax, m)`:
+Evaluate the associated Legendre Polynomials for order `m` and degree `l` in `lmin:lmax`
+as `collectPlm(x; m, lmin, lmax)`. By default `lmin` is chosen to be `abs(m)`, and may be omitted.
 
 ```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
+3-element OffsetArray(::Vector{Float64}, 3:5) with eltype Float64 with indices 3:5:
   -9.742785792574933
- -34.099750274012266
- -42.62468784251533
+ -34.09975027401223
+ -42.62468784251535
 ```
 
 Evaluate all the `n`th derivatives as `collectdnPl(x; lmax, n)`:
@@ -123,6 +155,8 @@ 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