FastTransforms submodule for the computation of some elliptic integrals and functions.

Complete elliptic integrals of the first and second kinds:

\[K(k) = \int_0^{\frac{\pi}{2}} \frac{{\rm d}\theta}{\sqrt{1-k^2\sin^2\theta}},\quad{\rm and},\]

\[E(k) = \int_0^{\frac{\pi}{2}} \sqrt{1-k^2\sin^2\theta} {\rm\,d}\theta.\]

Jacobian elliptic functions:

\[x = \int_0^{\operatorname{sn}(x,k)} \frac{{\rm d}t}{\sqrt{(1-t^2)(1-k^2t^2)}},\]

\[x = \int_{\operatorname{cn}(x,k)}^1 \frac{{\rm d}t}{\sqrt{(1-t^2)[1-k^2(1-t^2)]}},\]

\[x = \int_{\operatorname{dn}(x,k)}^1 \frac{{\rm d}t}{\sqrt{(1-t^2)(t^2-1+k^2)}},\]

and the remaining nine are defined by:

\[\operatorname{pq}(x,k) = \frac{\operatorname{pr}(x,k)}{\operatorname{qr}(x,k)} = \frac{1}{\operatorname{qp}(x,k)}.\]

Pre-plan an Inverse Padua Transform.

Pre-computes a 2D nonuniform fast Fourier transform.

For best performance, choose the right number of threads by FFTW.set_num_threads(4), for example.

Pre-computes a nonuniform fast Fourier transform of type N.

For best performance, choose the right number of threads by FFTW.set_num_threads(4), for example.

Pre-plan a Padua Transform.

Represent a scaled Toeplitz∘Hankel matrix:

DL(T∘H)DR

where the Hankel matrix H is non-negative definite, via

∑_{k=1}^r Diagonal(L[:,k])*T*Diagonal(R[:,k])

where L and R are determined by doing a rank-r pivoted Cholesky decomposition of H, which in low rank form is

H ≈ ∑_{k=1}^r C[:,k]C[:,k]'

so that L[:,k] = DL*C[:,k] and R[:,k] = DR*C[:,k].

This allows a Cholesky decomposition in 𝒪(K²N) operations and 𝒪(KN) storage, K = log N log ɛ⁻¹. The tuple storage allows plans applied to each dimension.

Pre-computes an inverse nonuniform fast Fourier transform of type N.

For best performance, choose the right number of threads by FFTW.set_num_threads(4), for example.

mpfr_t <: AbstractFloat

A Julia struct that exactly matches mpfr_t.

nearestjacobi_par(α, γ)

returns a number that is an integer different than γ but less than 1 away from α.

associatedjac2jac(v::AbstractVector, c::Integer, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false)

Convert the vector of expansions coefficients v from an associated Jacobi basis of orders (α,β) to a Jacobi basis of order (γ,δ). The keyword arguments denote whether the bases are normalized.

cheb2jac(v::AbstractVector, α, β; normcheb::Bool=false, normjac::Bool=false)

Convert the vector of expansions coefficients v from a Chebyshev basis to a Jacobi basis of order (α,β). The keyword arguments denote whether the bases are normalized.

cheb2leg(v::AbstractVector; normcheb::Bool=false, normleg::Bool=false)

Convert the vector of expansions coefficients v from a Chebyshev to a Legendre basis. The keyword arguments denote whether the bases are normalized.

cheb2ultra(v::AbstractVector, λ; normcheb::Bool=false, normultra::Bool=false)

Convert the vector of expansions coefficients v from a Chebyshev basis to an Ultraspherical basis of order λ. The keyword arguments denote whether the bases are normalized.

Modified Chebyshev moments of the first kind with respect to the absolute value weight:

\[ \int_{-1}^{+1} T_n(x) |x|{\rm\,d}x.\]

Modified Chebyshev moments of the second kind with respect to the absolute value weight:

\[ \int_{-1}^{+1} U_n(x) |x|{\rm\,d}x.\]

Modified Chebyshev moments of the first kind with respect to the Jacobi weight:

\[ \int_{-1}^{+1} T_n(x) (1-x)^\alpha(1+x)^\beta{\rm\,d}x.\]

Modified Chebyshev moments of the second kind with respect to the Jacobi weight:

\[ \int_{-1}^{+1} U_n(x) (1-x)^\alpha(1+x)^\beta{\rm\,d}x.\]

Modified Chebyshev moments of the first kind with respect to the logarithmic weight:

\[ \int_{-1}^{+1} T_n(x) \log\left(\frac{1-x}{2}\right){\rm\,d}x.\]

Modified Chebyshev moments of the second kind with respect to the logarithmic weight:

\[ \int_{-1}^{+1} U_n(x) \log\left(\frac{1-x}{2}\right){\rm\,d}x.\]

Modified Chebyshev moments of the first kind:

\[ \int_{-1}^{+1} T_n(x) {\rm\,d}x.\]

Modified Chebyshev moments of the second kind:

\[ \int_{-1}^{+1} U_n(x) {\rm\,d}x.\]

chebyshevtransform!(x, kind=Val(1))

transforms from values on a Chebyshev grid of the first or second kind to Chebyshev coefficients, in-place

chebyshevtransform(x, kind=Val(1))

transforms from values on a Chebyshev grid of the first or second kind to Chebyshev coefficients.

chebyshevutransform(x, ::Val{kind}=Val(1))

transforms from values on a Chebyshev grid of the first or second kind to Chebyshev coefficients of the 2nd kind (Chebyshev U expansion).

clenshaw!(c, A, B, C, x)

evaluates the orthogonal polynomial expansion with coefficients c at points x, where A, B, and C are AbstractVectors containing the recurrence coefficients as defined in DLMF, overwriting x with the results.

If c is a matrix this treats each column as a separate vector of coefficients, returning a vector if x is a number and a matrix if x is a vector.

clenshaw!(c, A, B, C, x, ϕ₀, f)

evaluates the orthogonal polynomial expansion with coefficients c at points x, where A, B, and C are AbstractVectors containing the recurrence coefficients as defined in DLMF and ϕ₀ is the zeroth polynomial, overwriting f with the results.

clenshaw!(c, x, f)

evaluates the first-kind Chebyshev (T) expansion with coefficients c at points x, overwriting f with the results.

clenshaw!(c, x)

evaluates the first-kind Chebyshev (T) expansion with coefficients c at points x, overwriting x with the results.

clenshaw(c, x)

evaluates the first-kind Chebyshev (T) expansion with coefficients c at the points x. x may also be a single Number.

Compute nodes of the Clenshaw—Curtis quadrature rule.

Compute weights of the Clenshaw—Curtis quadrature rule with modified Chebyshev moments of the first kind $\mu$.

Compute nodes of Fejer's first quadrature rule.

Compute nodes of Fejer's second quadrature rule.

Compute weights of Fejer's first quadrature rule with modified Chebyshev moments of the first kind $\mu$.

Compute weights of Fejer's second quadrature rule with modified Chebyshev moments of the second kind $\mu$.

forwardrecurrence!(v, A, B, C, x)

evaluates the orthogonal polynomials at points x, where A, B, and C are AbstractVectors containing the recurrence coefficients as defined in DLMF, overwriting v with the results.

Calculates the Gaunt coefficients in 64-bit floating-point arithmetic.

Calculates the Gaunt coefficients, defined by:

\[a(m,n,\mu,\nu,q) = \frac{2(n+\nu-2q)+1}{2} \frac{(n+\nu-2q-m-\mu)!}{(n+\nu-2q+m+\mu)!} \int_{-1}^{+1} P_n^m(x) P_\nu^\mu(x) P_{n+\nu-2q}^{m+\mu}(x) {\rm\,d}x.\]

or defined by:

\[P_n^m(x) P_\nu^\mu(x) = \sum_{q=0}^{q_{\rm max}} a(m,n,\mu,\nu,q) P_{n+\nu-2q}^{m+\mu}(x)\]

This is a Julia implementation of the stable recurrence described in:

Y.-l. Xu, Fast evaluation of Gaunt coefficients: recursive approach, J. Comp. Appl. Math., 85:53–65, 1997.

Compute a typed 0.5.

Computes an inverse nonuniform fast Fourier transform of type I.

Computes an inverse nonuniform fast Fourier transform of type II.

Inverse Padua Transform maps the 2D Chebyshev coefficients to the values of the interpolation polynomial at the Padua points.

jac2cheb(v::AbstractVector, α, β; normjac::Bool=false, normcheb::Bool=false)

Convert the vector of expansions coefficients v from a Jacobi basis of order (α,β) to a Chebyshev basis. The keyword arguments denote whether the bases are normalized.

jac2jac(v::AbstractVector, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false)

Convert the vector of expansions coefficients v from a Jacobi basis of order (α,β) to a Jacobi basis of order (γ,δ). The keyword arguments denote whether the bases are normalized.

jac2ultra(v::AbstractVector, α, β, λ; normjac::Bool=false, normultra::Bool=false)

Convert the vector of expansions coefficients v from a Jacobi basis of order (α,β) to an Ultraspherical basis of order λ. The keyword arguments denote whether the bases are normalized.

lag2lag(v::AbstractVector, α, β; norm1::Bool=false, norm2::Bool=false)

Convert the vector of expansions coefficients v from a Laguerre basis of order α to a La basis of order β. The keyword arguments denote whether the bases are normalized.

The principal branch of the Lambert-W function, defined by $x = W_0(x) e^{W_0(x)}$, computed using Halley's method for $x \in [-e^{-1},\infty)$.

leg2cheb(v::AbstractVector; normleg::Bool=false, normcheb::Bool=false)

Convert the vector of expansions coefficients v from a Legendre to a Chebyshev basis. The keyword arguments denote whether the bases are normalized.

maybereal(::Type{T}, x)

Return real-valued part of x if T is a type of a real number, and x otherwise.

modifiedherm2herm(v::AbstractVector{T}, u::Vector{T}; verbose::Bool=false)
modifiedherm2herm(v::AbstractVector{T}, u::Vector{T}, v::Vector{T}; verbose::Bool=false) where {T}

modifiedjac2jac(v::AbstractVector{T}, α, β, u::Vector{T}; verbose::Bool=false) where {T}
modifiedjac2jac(v::AbstractVector{T}, α, β, u::Vector{T}, v::Vector{T}; verbose::Bool=false) where {T}

modifiedlag2lag(v::AbstractVector{T}, α, u::Vector{T}; verbose::Bool=false)
modifiedlag2lag(v::AbstractVector{T}, α, u::Vector{T}, v::Vector{T}; verbose::Bool=false) where {T}

Computes a nonuniform fast Fourier transform of type I:

\[f_j = \sum_{k=0}^{N-1} c_k e^{-2\pi{\rm i} \frac{j}{N} \omega_k},\quad{\rm for}\quad 0 \le j \le N-1.\]

Computes a 2D nonuniform fast Fourier transform of type I-I:

\[F_{i,j} = \sum_{k=0}^{M-1}\sum_{\ell=0}^{N-1} C_{k,\ell} e^{-2\pi{\rm i} (\frac{i}{M} \omega_k + \frac{j}{N} \pi_{\ell})},\quad{\rm for}\quad 0 \le i \le M-1,\quad 0 \le j \le N-1.\]

Computes a nonuniform fast Fourier transform of type II:

\[f_j = \sum_{k=0}^{N-1} c_k e^{-2\pi{\rm i} x_j k},\quad{\rm for}\quad 0 \le j \le N-1.\]

Computes a 2D nonuniform fast Fourier transform of type II-II:

\[F_{i,j} = \sum_{k=0}^{M-1}\sum_{\ell=0}^{N-1} C_{k,\ell} e^{-2\pi{\rm i} (x_i k + y_j \ell)},\quad{\rm for}\quad 0 \le i \le M-1,\quad 0 \le j \le N-1.\]

Computes a nonuniform fast Fourier transform of type III:

\[f_j = \sum_{k=0}^{N-1} c_k e^{-2\pi{\rm i} x_j \omega_k},\quad{\rm for}\quad 0 \le j \le N-1.\]

Returns coordinates of the $(n+1)(n+2)/2$ Padua points.

Padua Transform maps from interpolant values at the Padua points to the 2D Chebyshev coefficients.

Creates $(n+2)x(n+1)$ matrix of interpolant values on the tensor grid at the $(n+1)(n+2)/2$ Padua points.

Vectorizes the function values at the Padua points.

Pre-computes an inverse nonuniform fast Fourier transform of type I.

Pre-computes an inverse nonuniform fast Fourier transform of type II.

Pre-plan an Inverse Padua Transform.

Pre-computes a 2D nonuniform fast Fourier transform of type I-I.

Pre-computes a nonuniform fast Fourier transform of type I.

Pre-computes a 2D nonuniform fast Fourier transform of type II-II.

Pre-computes a nonuniform fast Fourier transform of type II.

Pre-computes a nonuniform fast Fourier transform of type III.

Pre-plan a Padua Transform.

Pochhammer symbol $(x)_n = \frac{\Gamma(x+n)}{\Gamma(x)}$ for the rising factorial.

BigFloat is a mutable struct and there is no guarantee that each entry in an AbstractArray{BigFloat} is unique. For example, looking at the Limbs,

Id = Matrix{BigFloat}(I, 3, 3)
map(x->x.d, Id)

shows that the ones and the zeros all share the same pointers. If a C function assumes unicity of each datum, then the array must be renewed with a deepcopy.

Pointwise evaluation of real orthonormal spherical harmonic:

\[Y_\ell^m(\theta,\varphi) = (-1)^{|m|}\sqrt{(\ell+\frac{1}{2})\frac{(\ell-|m|)!}{(\ell+|m|)!}} P_\ell^{|m|}(\cos\theta) \sqrt{\frac{2-\delta_{m,0}}{2\pi}} \left\{\begin{array}{ccc} \cos m\varphi & {\rm for} & m \ge 0,\\ \sin(-m\varphi) & {\rm for} & m < 0.\end{array}\right.\]

Stirling's asymptotic series for $\Gamma(z)$.

Creates $(n+2)x(n+1)$ Chebyshev coefficient matrix from triangle coefficients.

Creates length $(n+1)(n+2)/2$ vector from matrix of triangle Chebyshev coefficients.

Pointwise evaluation of triangular harmonic:

\[\tilde{P}_{\ell,m}^{(\alpha,\beta,\gamma)}(x,y).\]

Compute a typed 2.

ultra2cheb(v::AbstractVector, λ; normultra::Bool=false, normcheb::Bool=false)

Convert the vector of expansions coefficients v from an Ultraspherical basis of order λ to a Chebyshev basis. The keyword arguments denote whether the bases are normalized.

ultra2jac(v::AbstractVector, λ, α, β; normultra::Bool=false, normjac::Bool=false)

Convert the vector of expansions coefficients v from an Ultraspherical basis of order λ to a Jacobi basis of order (α,β). The keyword arguments denote whether the bases are normalized.

ultra2ultra(v::AbstractVector, λ, μ; norm1::Bool=false, norm2::Bool=false)

Convert the vector of expansions coefficients v from an Ultraspherical basis of order λ to an Ultraspherical basis of order μ. The keyword arguments denote whether the bases are normalized.

For 64-bit floating-point arithmetic, the Lambda function uses the asymptotic series for $\tau$ in Appendix B of

I. Bogaert and B. Michiels and J. Fostier, 𝒪(1) computation of Legendre polynomials and Gauss–Legendre nodes and weights for parallel computing, SIAM J. Sci. Comput., 34:C83–C101, 2012.

The Lambda function $\Lambda(z) = \frac{\Gamma(z+\frac{1}{2})}{\Gamma(z+1)}$ for the ratio of gamma functions.

The Lambda function $\Lambda(z,λ₁,λ₂) = \frac{\Gamma(z+\lambda_1)}{Γ(z+\lambda_2)}$ for the ratio of gamma functions.

The Kronecker $\delta$ function:

\[\delta_{k,j} = \left\{\begin{array}{ccc} 1 & {\rm for} & k = j,\\ 0 & {\rm for} & k \ne j.\end{array}\right.\]