| Up | Next | Prev | PrevTail | Tail |
This special function package is separated into two portions to make it easier to handle. The packages are called SPECFN and SPECFN2. The first one is more general in nature, whereas the second is devoted to special special functions. Additional examples can be found in the files specfn.tst and specfn2.tst in the packages/specfn directory.
Authors: Chris Cannam, with contributions from Winfried Neun, Herbert Melenk, Victor Adamchik, Francis Wright, Alan Barnes and several others.
\(\newcommand {\f }[1]{\texttt {#1}}\)The package SPECFN is designed to provide algebraic and numeric manipulations of many common special functions, namely:
The psi function & its derivatives;
The Bessel functions \(J\) and \(Y\) of the first and second kinds;
The modified Bessel functions \(I\) and \(K\);
The Hankel functions \(H^{(1)}\) and \(H^{(2)}\);
The Kummer hypergeometric functions M and U;
Associated Legendre Functions (Spherical and Solid Harmonics);
and some well-known constants.
All of the above functions (except Stirling numbers) are autoloading.
More information on all these functions may be found on the website DLMF:NIST although currently not all functions may conform to these standards.
All algorithms whose sources are uncredited are culled from series or expressions found in the Dover Handbook of Mathematical Functions[AS72].
A nice collection of example plots for special functions can be found in the file specplot.tst in the subfolder plot of the packages folder. These examples will reproduce a number of well-known pictures from [AS72].
Most of these polynomial functions are not autoloading. This package needs to be loaded before they may be used with the command:
load_package specfn;
The polynomial function sets available are:
All of the operators supported by this package have certain algebraic simplification rules to handle special cases, poles, derivatives and so on. Such rules are applied whenever they are appropriate. However, if the rounded switch is on, numeric evaluation is also carried out. Unless otherwise stated below, the result of an application of a special function operator to real or complex numeric arguments in rounded mode will be approximated numerically whenever it is possible to do so. All approximations are to the current precision.
Most algebraic simplifications within the special function package are defined in the form of a REDUCE ruleset. Therefore, in order to get a quick insight into the simplification rules one can use the ShowRules operator, e.g.
ShowRules BesselI;
1 ~z - ~z
{besseli(~n,~z) => ---------------*(e - e )
sqrt(pi*2*~z)
1
when numberp(~n) and ~n=---,
2
1 ~z - ~z
besseli(~n,~z) => ---------------*(e + e )
sqrt(pi*2*~z)
1
when numberp(~n) and ~n= - ---,
2
besseli(~n,~z) => 0
when numberp(~z) and ~z=0
and numberp(~n) and ~n neq 0,
besseli(~n,~z) => besseli( - ~n,~z) when numberp(~n)
and impart(~n)=0 and ~n=floor(~n) and ~n<0,
besseli(~n,~z) => do*i(~n,~z)
when numberp(~n) and numberp(~z) and *rounded,
df(besseli(~n,~z),~z)
besseli(~n - 1,~z) + besseli(~n + 1,~z)
=> -----------------------------------------,
2
df(besseli(~n,~z),~z)
=> besseli(1,~z) when numberp(~n) and ~n=0}
Several REDUCE packages (such as Sum or Limits) obtain different (hopefully better) results for the algebraic simplifications when the SPECFN package is loaded, because the latter package contains some information which may be useful and directly applicable for other packages, e.g.:
sum(1/k^s,k,1,infinity); % evaluates to zeta(s)
A record is kept of all values previously approximated, so that should a value be required which has already been computed to the current precision or greater, it can be simply looked up. This can result in some storage overheads, particularly if many values are computed which will not be needed again. In this case, the switch savesfs may be turned off in order to inhibit the storage of approximated values. The switch is on by default.
The SPECFN package includes manipulation and limited numerical evaluation for some integral functions, namely
erf, erfc, Si, Shi, si, Ci, Chi, Ei, Li, Fresnel_C, and Fresnel_S.
The error function, its complement and the two Fresnel integrals are defined by: \begin {align*} \mathrm {erf}(z) & = \frac {2}{\sqrt \pi }\int _0^z e^{-t^2}\, \mathrm {d}t\\ \mathrm {erfc}(z) & = \frac {2}{\sqrt \pi }\int _z^\infty e^{-t^2}\, \mathrm {d}t = 1 - \mathrm {erf}(z)\\ C(z) & = \int _0^z \cos \left (\frac {\pi }{2} t^2\right )\, \mathrm {d}t\\ S(z) & = \int _0^z \sin \left (\frac {\pi }{2} t^2\right )\, \mathrm {d}t \end {align*}
respectively.
The exponential and related integrals are defined by the following: \begin {align*} \mathrm {Ei}(z) & = e^{-z}\int _z^\infty \frac {e^{-t}} {t+z}\, \mathrm {d}t\\ \mathrm {Li}(z) & = \int _0^z \frac {\mathrm {d}t} {\log t}\\ \mathrm {Si}(z) & = \int _0^z \frac {\sin t} {t}\, \mathrm {d}t\\ \mathrm {si}(z) & = -\int _z^\infty \frac {\sin t} {t}\, \mathrm {d}t = \mathrm {Si}(z) - \frac {\pi }{2}\\ \mathrm {Ci}(z) & = -\int _z^\infty \frac {\cos t} {t}\, \mathrm {d}t = \int _0^z \frac {\cos t -1} {t}\, \mathrm {d}t + \log z +\gamma \\ \mathrm {Shi}(z) & = \int _0^z \frac {\sinh t} {t}\, \mathrm {d}t\\ \mathrm {Chi}(z) & = \int _0^z \frac {\cosh t -1} {t}\, \mathrm {d}t +\log z +\gamma \end {align*}
where \(\gamma \) is Euler’s constant (Euler_gamma).
The definitions of the exponential and related integrals, the derviatives and some limits are known, together with some simple properties such as symmetry conditions.
The numerical approximations for the integral functions suffer from the fact that the precision is not set correctly for values of the argument above 10.0 (approx.) and from the usage of summations even for large arguments.
\(\mathop {\mathrm {Li}(z)}\) is simplified to \(\mathop {\mathrm {Ei}}(\ln (z))\) .
This is represented by the unary operator Gamma. The Gamma function is defined by the integral: \[ \Gamma (a) = \int _0^\infty e^{-t}t^{a-1}\, \mathrm {d}t.\] Initial transformations applied with rounded off are: \(\Gamma (n)\) for integral \(n\) is computed, \(\Gamma (n+1/2)\) for integral \(n\) is rewritten to an expression in \(\sqrt \pi \), \(\Gamma (n+1/m)\) for natural \(n\) and \(m\) a positive integral power of 2 less than or equal to 64 is rewritten to an expression in \(\Gamma (1/m)\), expressions with arguments at which there is a pole are replaced by infinity, and those with a negative (real) argument are rewritten so as to have positive arguments.
The algorithm used for numerical approximation is an implementation of an asymptotic series for \(\ln (\Gamma )\), with a scaling factor obtained from the Pochhammer symbols.
An expression for \(\Gamma '(z)\) in terms of \(\Gamma \) and \(\psi \) is included.
The (unnormalised) incomplete gamma function is provided by the binary function m_gamma. In the literature it is normally represented as \(\gamma (a,z)\) and is defined by \[ \gamma (a, z) = \int _0^z e^{-t}t^{a-1}\, \mathrm {d}t.\] The normalised incomplete gamma function \(P(a,z)\) is provided by the binary function igamma and is defined as \[P(a,z) = \frac {\gamma (a,z)}{\Gamma (a)}.\]
The binary function \(B(a,b)\) is related to the \(\Gamma \) function[AS72] and is defined by \[ B(a,b) = \int _0^1 t^a (1-t)^b\, \mathrm {d}t = \frac {\Gamma (a)\Gamma (b)}{\Gamma (a+b)}.\] It is represented by the binary function Beta.
The unnormalised and nomalised incomplete Beta funtions are defined by \begin {align*} B_x(a,b) & = \int _0^x t^a (1-t)^b\, \mathrm {d}t,\\ I_x(a,b) & = \frac {B_x(a,b)}{B(a,b)} \end {align*}
respectively. \(I_x(a,b)\) is represented by the ternary function ibeta(a,b,x).
This is represented by the unary operator psi. It is defined as the logarithmic derivative of the \(\Gamma \) function: \[ \psi (z) = \frac {\Gamma '(z)}{\Gamma (z)}. \]
Initial transformations for \(\psi \) are applied on a similar basis to those for \(\Gamma \); where possible, \(\psi (x)\) is rewritten in terms of \(\psi (1)\) and \(\psi (\frac {1}{2})\), and expressions with negative arguments are rewritten to have positive ones.
The algorithm for numerical evaluation of \(\psi \) is based upon an asymptotic series, with a suitable scaling.
Relations for the derivative and integral of \(\psi \) are included.
The \(n\)th derivative of the \(\psi \) function is represented by the binary operator Polygamma, whose first argument is \(n\).
Initial manipulations on \(\psi ^{(n)}\) are few; where the second argument is \(1\) or \(3/2\), the expression is rewritten to one involving the Riemann \(\zeta \) function, and when the first is zero it is rewritten to \(\psi \); poles are also handled.
Numerical evaluation is available for real and complex arguments. The algorithm used is again an asymptotic series with a scaling factor; for negative (second) arguments, a Reflection Formula is used, introducing a term in the \(n\)th derivative of \(\cot (z\pi )\).
Simple relations for derivatives and integrals are provided.
Support is provided for the Bessel functions \(J\) and \(Y\), the modified Bessel functions \(I\) and \(K\), and the Hankel functions of the first and second kinds. The relevant operators are, respectively, BesselJ, BesselY, BesselI, BesselK, Hankel1 and Hankel2, which are all binary operators.
The Bessel functions \(J_\nu (z)\) and \(Y_\nu (z)\) are solutions of the Bessel equation: \[z^2\frac {\mathrm {d}^2w}{\mathrm {d}z^2}+z\frac {\mathrm {d}w}{\mathrm {d}z} + (z^2 - \nu ^2)w = 0.\] Bessel’s function of the first kind, \(J_\nu (z)\), has the series expansion: \[J_\nu (z) = \left (\frac {z}{2}\right )^\nu \sum _{k=0}^\infty (-1)^k \frac {(z/2)^{2k}}{k! \Gamma (\nu + k+1)}\,.\] Bessel’s function of the second kind, \(Y_\nu (z)\), (for non-integral \(\nu \)) is defined by: \[Y_\nu (z) = \frac {J_\nu (z) \cos (\nu \pi )-J_{-\nu }(z)}{\sin (\nu \pi )}\] or by its limiting value: \[Y_\nu (z) = \left .\frac {1}{\pi }\frac {\partial {J_\nu (z)}}{\partial \nu }\right |_{\nu =n} +\quad \left .\frac {(-1)^n}{\pi }\frac {\partial {J_\nu (z)}}{\partial \nu }\right |_{\nu =-n}. \] It is sometimes known as Weber’s function.
The Hankel functions are alternative solutions of the Bessel equation distinguished by their asymptotic behaviour as \(z\rightarrow \infty \): \begin {align*} H_\nu ^{(1)}(z) & \sim \sqrt {\frac {2}{\pi z}}\exp \left (i\left (z-\frac {\nu \pi } {2}-\frac {\pi }{4}\right )\right ),\\ H_\nu ^{(2)}(z) & \sim \sqrt {\frac {2}{\pi z}}\exp \left (-i\left (z-\frac {\nu \pi } {2}-\frac {\pi }{4}\right )\right ). \end {align*}
The modified Bessel functions \(I_\nu (z)\) and \(K_\nu (z)\) are solutions of the modified Bessel equation: \[z^2\frac {\mathrm {d}^2w}{\mathrm {d}z^2}+z\frac {\mathrm {d}w}{\mathrm {d}z} - (z^2 + \nu ^2)w = 0\,.\] Since they may be obtained by replacing \(z\) by \(\pm i z\) the modified Bessel functions are sometimes called Bessel functions of imaginary argument. \(I_\nu (z)\) has the series expansion: \[I_\nu (z) = \left (\frac {z}{2}\right )^\nu \sum _{k=0}^\infty \frac {(z/2)^{2k}}{k!\Gamma (\nu +k+1)}\,,\] whereas \(K_\nu (z)\) is distinguished by its asymptotic behaviour: \[K_\nu (z) \sim \sqrt {\frac {\pi }{2 z}}e^{-z}\] as \(z\rightarrow \infty \). For more information, see the DLMF:NIST chapters on Hankel & Bessel functions and Modified Bessel functions.
The following initial transformations are performed:
trivial cases or poles of \(J\), \(Y\), \(I\) and \(K\) are handled;
\(J\), \(Y\), \(I\) and \(K\) with negative first argument are transformed to have positive first argument;
\(J\) with negative second argument is transformed to have positive second argument;
\(Y\) or \(K\) with non-integral or complex second argument is transformed into an expression in \(J\) or \(I\) respectively;
derivatives of \(J\), \(Y\) and \(I\) are carried out;
derivatives of \(K\) with zero first argument are carried out;
derivatives of Hankel functions are carried out.
Also, if the complex switch is on and rounded is off, expressions in Hankel functions are rewritten in terms of Bessel functions.
No numerical approximation is provided for the Bessel \(K\) function, or for the Hankel functions for anything other than special cases. The algorithms used for the other Bessel functions are generally implementations of standard ascending series for \(J\), \(Y\) and \(I\), together with asymptotic series for \(J\) and \(Y\); usually, the asymptotic series are tried first, and if the argument is too small for them to attain the current precision, the standard series are applied. An obvious optimization prevents an attempt with the asymptotic series if it is clear from the outset that it will fail.
There are no rules for the integration of Bessel and Hankel functions.
Support is provided for the Airy Functions \(Ai\) and \(Bi\) and for their derivatives \(Ai'\) and \(Bi'\). The relevant operators are respectively Airy_Ai, Airy_Bi, Airy_Aiprime and Airy_Biprime, which are all unary.
Airy functions are solutions of the differential equation: \[ \frac {\mathrm {d}^2 w}{\mathrm {d}z^2} = z w.\]
Trivial cases of Airy_Ai and Airy_Bi and their primes are evaluated, and all functions accept both real and complex arguments.
The Airy Functions can also be represented in terms of Bessel Functions by activating an inactive rule set:
let Airy2Bessel_rules;
As a result the Airy_Ai function will be evaluated using the formula: \[ \f {Ai}(z) = \frac {1}{3} \sqrt {z} \left [I_{-1/3}(\zeta ) - I_{1/3}(\zeta )\right ], \text { where } \zeta = \frac {2}{3} z^{\frac {2}{3}}. \] Note: In order to obtain satisfactory approximations to numerical values both the complex and rounded switches must be on.
The algorithms used for the Airy Functions are implementations of standard ascending series, together with asymptotic series. At some point it is better to use the asymptotic rather than the ascending series, which is calculated by the program and depends on the given precision.
There are no rules for the integration of Airy Functions.
This package also provides some support for other functions, in the form of algebraic simplifications:
The Struve H and L functions, through the binary operators StruveH and StruveL, for which manipulations are provided to handle special cases, simplify to more readily handled functions where appropriate, and differentiate with respect to the second argument. These functions with arguments \(\nu \) and \(x\) are solutions of the differential equation: \[\frac {\mathrm {d}^2w}{\mathrm {d}x^2}+\frac {1}{x}\frac {\mathrm {d}w}{\mathrm {d}x} +\left (1-\frac {\nu ^2}{x^2}\right )w = \frac {(z/2)^{\nu -1}}{\sqrt {\pi }\Gamma (\nu +1/2)}\,.\]
The Lommel functions of the first and second kinds, through the ternary operators Lommel1 and Lommel2 with arguments \(\nu \), \(\mu \) and \(x\) may be considered generalisations of the Struve functions satisfying the differential equation: \[\frac {\mathrm {d}^2w}{\mathrm {d}x^2}+\frac {1}{x}\frac {\mathrm {d}w}{\mathrm {d}x} +\left (1-\frac {\nu ^2}{x^2}\right )w = z^{\mu -1}\,.\] Manipulations are provided to handle special cases and simplify where appropriate.
The Kummer confluent hypergeometric functions M and U (the hypergeometric \(_1F_1\) or \(\Phi \), and \(z^{-a}{_2F_0}\) or \(\Psi \), respectively), represented by the ternary operators KummerM and KummerU with arguments \(a\), \(b\) and \(x\), are solutions of the differential equation: \[\frac {\mathrm {d}^2w}{\mathrm {d}x^2}+(b-x)\frac {\mathrm {d} w}{\mathrm {d}x} -a w = 0\,.\] There are manipulations for special cases and simplifications, derivatives and, for the M function, numerical approximations for real arguments.
The Whittaker M and W functions are variations upon the Kummer functions, which are represented by the ternary operators WhittakerM and WhittakerW with arguments \(\kappa \), \(\mu \) and \(x\). They satisfy the Whittaker differential equation: \[\frac {\mathrm {d}^2W}{\mathrm {d}x^2}+\left (\frac {1-4\mu ^2}{4x^2}+ \frac {\kappa }{x} -\frac {1}{4}\right )W =0\,,\] which is obtained from the Kummer differential equation via the substituions \[W=e^{z/2}z^{\mu +1/2} w,\qquad \kappa =b/2-a\qquad \mu =(b-1)/2\,.\] The Whittaker M and W functions with non-numeric arguments are simplified to expressions involving the Kummer M and U functions respectively.
This is represented by the unary operator Zeta and defined by the formula: \[\zeta (s) = \sum _{n=1}^\infty \frac {1}{n^s}.\]
With rounded off, \(\zeta (z)\) is evaluated numerically for even integral arguments in the range \(-31 < z < 31\), and for odd integral arguments in the range \(-30 < z < 16\). Outside this range the values become a little unwieldy.
Numerical evaluation of \(\zeta \) is only carried out if the argument is real. The algorithms used for \(\zeta \) are: for odd integral arguments, an expression relating \(\zeta (n)\) with \(\psi ^{n-1}(3)\); for even arguments, a trivial relationship with the Bernoulli numbers; and for other arguments the approach is either (for larger arguments) to take the first few primes in the standard over-all-primes expansion, and then continue with the defining series with natural numbers not divisible by these primes, or (for smaller arguments) to use a fast-converging series obtained from [BO78].
There are no rules for differentiation or integration of \(\zeta \).
The dilogarithm function \(Li_2(z)\) is defined by \[\mathrm {Li}_2(z) \equiv \sum _{n=1}^\infty \frac {z^n}{n^2} = -\int _0^z \frac {\log (1-t)}{t}\, \mathrm {d}t\] and represented by the unary function dilog.
The polylogarithm function \(Li_s(z)\) is defined by \[\mathrm {Li}_s(z) \equiv \sum _{n=1}^\infty \frac {z^n}{n^s} = \frac {z}{\Gamma (s)} \int _0^\infty \frac {t^{s-1}}{e^t-z}\, \mathrm {d}t.\] and represented by the binary function Polylog. The case \(s=2\) is, of course, the dilogarithm function and the special case when \(z=1\) gives the Riemann zeta function \(\zeta (s)\). For \(s=1\), the polylogarithm reduces to the elementary function: \(-\log (1-t)\).
Lerch’s transcendent or Lerch Phi function is defined by \[\Phi (z,s,a) = \sum _{n=0}^\infty \frac {z^n}{(n+a)^s}.\] It is represented by the ternary function Lerch_Phi(z,s,a). For the special case \(a=1\), Lerch’s function is related to a polylogarithm: \(z Li_s(z) = \Phi (z,s,1)\).
Lambert’s function \(\omega (x)\), represented by the unary operator Lambert_W, is the inverse of the function \(x=we^w\). Therefore it is an important contribution for the solve package.
For real-valued arguments \(\omega (x)\) is only real-valued in the interval \((-1/e, \infty )\). In the interval \((-1/e, 0)\), it is double-valued with a branch point at the point (-1/e, -1) where \(\omega '(x)\) is singular. The positive branch is defined on the interval \((-1/e, \infty )\) where it is monotonically increasing with \(\omega (x) > -1\). The negative branch is defined on the interval \((-1/e, 0)\) where it is monotonically decreasing with \(\omega (x) < -1\).
Simplification rules for \(\omega (x)\) are provided for the special arguments \(0\) and \(-1/e\) and for its logarithm, derivative and integral. A previous rule for its exponential caused problems with power series expansions about zero and has been deactivated. This does not seem to impact on the SOLVE package. However, this rule may be reactivated if required by
let lambert_exp_rule; % and deactivated again by clear lambert_exp_rule;
The function is studied extensively in [HCGJ92]. The current implementation will compute values on the principal branch for all complex numerical arguments only if the switch rounded is on. However, since the numerical computations are carried out in complex-rounded mode, it is also better to turn the switch complex to on to avoid repeated irritating mode change warnings.
The real positive branch is part of the principal branch and currently there is no way of computing values on the real negative branch or indeed any non-principal values.
The relevant operators are, respectively,
SolidHarmonicY and SphericalHarmonicY.
The SolidHarmonicY operator implements the Solid Harmonics described below. It expects 6 parameter, namely \(n\), \(m\), \(x\), \(y\), \(z\) and \(r2\) and returns a polynomial in \(x\), \(y\), \(z\) and \(r2\).
The operator SphericalHarmonicY is a special case of SolidHarmonicY with the usual definition:
algebraic procedure SphericalHarmonicY(n,m,theta,phi); SolidHarmonicY(n,m,sin(theta)*cos(phi), sin(theta)*sin(phi),cos(theta),1)$
Solid Harmonics of order \(n\) (Laplace polynomials) are homogeneous polynomials of degree \(n\) in \(x\), \(y\), \(z\) which are solutions of the Laplace equation:-
df(P,x,2) + df(P,y,2) + df(P,z,2) = 0.
There are \(2n+1\) independent such polynomials for any given \(n \geq 0\) and with:-
w!0 = z, w!+ = i*(x-i*y)/2, w!- = i*(x+i*y)/2,
they are given by the Fourier integral:-
S(n,m,w!-,w!0,w!+) = (1/(2*pi)) * for u:=-pi:pi integrate(w!0 + w!+ * exp(i*u) + w!- * exp(-i*u))^n * exp(i*m*u) * du;
which is obviously zero if \(|m| > n\) since then all terms in the expanded integrand contain the factor \(e^{iku}\) with \(k \neq 0\).
\(S(n,m,x,y,z)\) is proportional to
r^n * Legendre(n,m,cos theta) * exp(i*phi)
where \(r^2 = x^2 + y^2 + z^2\).
The spherical harmonics are simply the restriction of the solid harmonics to the surface of the unit sphere and the set of all spherical harmonics with \(n \geq 0, -n \leq m \leq n\) form a complete orthogonal basis on it, i.e. \(\langle n,m | n',m' \rangle = \delta _{n,n'} \delta _{m,m'}\) using \(\langle \ldots | \ldots \rangle \) to designate the scalar product of functions over the spherical surface.
The coefficients of the solid harmonics are normalised in what follows to yield an orthonormal system of spherical harmonics.
Given their polynomial nature, there are many recursions formulae for the solid harmonics and any recursion valid for Legendre functions can be ‘translated’ into solid harmonics. However the direct proof is usually far simpler using Laplace’s definition.
It is also clear that all differentiations of solid harmonics are trivial, qua polynomials.
Some substantial reduction in the symbolic form would occur if one maintained throughout the recursions the symbol \(r2\) (\(r\) cannot occur as it is not rational in \(x,y,z\)). Formally the solid harmonics appear in this guise as more compact polynomials in \(x,y,z,r2\).
Only two recursions are needed:-
(i) along the diagonal \((n,n)\);
(ii) along a line of constant \(n\): \((m,m),(m+1,m),\ldots ,(n,m)\).
Numerically these recursions are stable.
For \(m < 0\) one has:- \[ S(n,m,x,y,z) = (-1)^m S(n,-m,x,-y,z). \]
The operators ThreeJSymbol and Clebsch_Gordan are defined as in [LB68] or [Edm57] and expect as arguments three lists of values {\(j_i,m_i\)}, e.g.
ThreeJSymbol({J+1,M},{J,-M},{1,0});
Clebsch_Gordan({2,0},{2,0},{2,0});
The operator SixJSymbol is defined as in [LB68] or [Edm57] and expects two lists of values {\(j_1,j_2,j_3\)} and {\(l_1,l_2,l_3\)} as arguments, e.g.
SixJSymbol({7,6,3},{2,4,6});
In the current implementation of SixJSymbol there is only limited reasoning about the minima and maxima of the summation using the INEQ package, such that in most cases the special 6j-symbols (see e.g. [LB68]) will not be found.
The Stirling numbers of the first and second kind are computed by calling the binary operators Stirling1 and Stirling2 respectively.
Stirling numbers of the first kind have the generating function: \[\sum _{m=0}^n s_n^m x^m = (x-n+1)_n\] where \((x-n+1)_n\) is the Pochhammer symbol. This provides a convenient way of calculating these Stirling numbers by extracting coefficients of the polynomial obtained by evaluating the Pochhammer symbol. REDUCE however uses an explicit summation.
Stirling numbers of the second kind are defined by the formula: \[S_n^m = \frac {1}{m!} \sum _{k=0}^m (-1)^{m-k} \binom {m}{k} k^n.\] REDUCE uses this explicit summation to evaluate Stirling numbers of the second kind.
The following well-known constants are defined in the REDUCE core, but the code for computing their numerical value when the switch ROUNDED is on is contained in the special function package.
Euler_Gamma : Euler’s constant, also available as \(-\psi (1)\);
Catalan : Catalan’s constant;
Khinchin : Khinchin’s constant, defined in [Khi64] (which takes a lot of time to compute);
Golden_Ratio : \(\displaystyle \frac {1 + \sqrt {5}}{2}\)
All the polynomials in this section take two or more parameters; the first is the degree of the polynomial and the last is its argument. Any remaining arguments are parameters which in the literature are normally rendered as subscripts and superscripts. First, the definitions appropriate to all the sets of orthogonal polynomials in the following subsections are listed.
A set of polynomials \(\{p_n(x)\},\, n=0,1,\ldots \) are said to be orthogonal on open interval \((a, b)\) (where \(a\) and/or \(b\) may be infinite) with positive weight function \(w(x)\) if \[\int _a^b p_n(x)p_m(x)w(x) \mathrm {d}x = 0 \qquad \mbox {when}\ m \neq n.\] This defines each polynomial \(p_n(x)\) up to a constant factor \(c_n\) which is usually fixed by normalisation. If these factors are chosen so that \[h_n=\int _a^b (p_n(x))^2 w(x) \mathrm {d}x = 1\qquad \mbox {i.e.}\ c_n=\sqrt {h_n}\] then the polynomial set is said to be orthonormal. An alternative normalisation, that is sometimes used, is to set the leading term of each polynomial \(k_n = 1\). The polynomial set is then said to be monic.
In REDUCE the normalisation is chosen so that the polynomial sets are orthonormal and hence \(k_n \neq 1\) in general. In the subsections below on each of the polynomial sets, the interval \((a,b)\) over which the polynomials are orthogonal, the weight function \(w(x)\) and the leading coefficient \(k_n\) of the polynomial of degree \(n\) are given together with any constraints on any additional parameters. Also given are what might be called the ‘first moment’ \(\tilde {h}_n\) of the \(n\)th polynomial defined by: \[\tilde {h}_n = \int _a^b x (p_n(x))^2 w(x) \,\mathrm {d}x\] and the ratio \[r_n = \frac {\tilde {k}_n}{k_n}\quad \mbox {where}\ p_n(x) = k_n x^n + \tilde {k}_n x^{n-1} \ldots \] These quantities may be used in recurrence relations when generating the polynomials.
The function call LegendreP(n,x) will return the \(n\)th Legendre polynomial if \(n\) is a non-negative integer; otherwise the result will involve the original operator LegendreP or on graphical interfaces \(P_n(x)\) will be output.
The interval of definition is \((-1, 1)\), the weight function \(w(x)=1\) and, for the orthonormal case, the leading coefficients are given by \(k_n=2^n (\frac {1}{2})_n/n!\) where \((\frac {1}{2})_n\) is the Pochhammer symbol. Also \(\tilde {h}_n = \frac {2}{2 n +1}\) and \(r_n =0\).
The function call LegendreP(n,m,x) will return the \(n\)th associated Legendre function if \(n\) and \(m\) are integers with \(0 \leq m \leq n\); otherwise the result will involve the original operator LegendreP or on graphical interfaces \(P_n^{(m)}(x)\) will be output.
They are defined by \[P_n^{(m)}(x) = (-1)^m(1-x^2)^{m/2}\frac {\mathrm {d}^m P_n(x)}{\mathrm {d}x^m}\,;\] it should be noted that they are only polynomials if \(m\) is even. Currently the extension of these functions to negative \(n\) and \(m\) is not implemented in REDUCE.
For fixed \(m\) these functions are orthogonal over the interval \((-1, 1)\); the weight function being \(w(x)=1\). However, unlike the polynomials in the rest of this section, they are not orthonormal: \[\int _{-1}^1 \left (P_n^{(m)}(x)\right )^2 \mathrm {d}x = h_n = \frac {2(l+m)!}{(2l+1)!(l-m)!}\,.\]
The function call ChebyshevT(n,x) will return the \(n\)th Chebyshev polynomial of the first kind if \(n\) is a non-negative integer; otherwise the result will involve the original operator ChebyshevT or on graphical interfaces \(T_n(x)\) will be output.
The interval of definition is \((-1, 1)\), the weight function \(w(x)=(1-x^2)^{-1/2}\) and, for the orthonormal case, the leading coefficients are given by \(k_n= 2^{n-1}\ \mbox {for}\ n>0;\ k_0 =1\). Also \(\tilde {h}_n = \pi /2\ \ \mbox {for}\ n>0;\ \tilde {h}_0 =\pi \,\) and \(r_n=0\).
The function call ChebyshevU(n,x) will return the \(n\)th Chebyshev polynomial of the second kind if \(n\) is a non-negative integer; otherwise the result will involve the original operator ChebyshevU or on graphical interfaces \(U_n(x)\) will be output.
The interval of definition is \((-1, 1)\), the weight function \(w(x)=(1-x^2)^{-1/2}\) and, for the orthonormal case, the leading coefficients are given by \(k_n= 2^n\), \(\tilde {h}_n = \pi /2\) and \(r_n=0\).
The function call GegenbauerP(n,a,x) will return the Gegenbauer polynomial of degree \(n\) and parameter \(a\) if \(n\) is a non-negative integer and \(a\) is numerical; otherwise the result will involve the original operator GegenbauerP or on graphical interfaces \(C_n^{(a)}(x)\) will be output.
The interval of definition is \((-1, 1)\), the weight function \(w(x)=(1-x^2)^{a-1/2}\) and, for the orthonormal case, the leading coefficients are given by \(k_n= 2^n (a)_n/n!\) where \((a)_n\) is the Pochhammer symbol. The parameter \(a\) should satisfy \(a >-1/2,\ \ a \neq 0\). Also \[\tilde {h}_n = \frac {2^{1-2 a}\pi \Gamma (n+2a)}{(n+a)(\Gamma (a))^2n!}\qquad \mbox {and}\quad r_n=0\,.\]
The function call JacobiP(n,a,b,x) will return the Jacobi polynomial of degree \(n\) and parameters \(a\) and \(b\) if \(n\) is a non-negative integer and \(a\) and \(b\) are numerical; otherwise the result will involve the original operator JacobiP or on graphical interfaces \(P_n^{(a, b)}(x)\) will be output.
The interval of definition is \((-1, 1)\), the weight function \(w(x)=(1-x)^a(1+x)^b\) and, for the orthonormal case, the leading coefficients are given by \[\tilde {h}_n = \frac {(n+a+b+1)_n}{2^nn!}\] where \((n+a+b+1)_n\) is the Pochhammer symbol. The parameters \(a\) and \(b\) should satisfy \(a >-1,\ \ b > -1\). Also \begin {align*} \tilde {h}_0 & = 2^{a+b+1}\frac {\Gamma (a+1)\Gamma (b+1)}{\Gamma (a+b+2)}\\ \tilde {h}_n & = 2^{a+b+1}\frac {\Gamma (n+a+1)\Gamma (n+b+1)}{(2 n +a+b+1) \Gamma (n+a+b+1)n!}\quad \mbox {for}\ n>0\\ r_n & = \frac {n(a-b)}{2n+a+b}. \end {align*}
The Legendre, Chebyshev and Gegenbauer polynomials are all, in fact, special cases of the Jacobi polynomials.
The function call LaguerreP(n,x) will return the \(n\)th Laguerre polynomial if \(n\) is a non-negative integer; otherwise the result will involve the original operator LaguerreP or on graphical interfaces \(L_n(x)\) will be output.
The interval of definition is \((0, \infty )\), the weight function \(w(x)=e^{-x}\) and, for the orthonormal case, the leading coefficients are given by \(k_n=(-1)^n/n!\), \(\tilde {h}_n =1\) and \(r_n = -n^2\).
If used with three arguments LaguerreP(n,a,x) returns the \(n\)th generalised (or associated) Laguerre polynomial if \(n\) is a non-negative integer and \(a\) is numeric; otherwise the result will involve the original operator LaguerreP or on graphical interfaces \(L_n^{(a)}(x)\) will be output. These are more properly called Sonin polynomials after their discoverer N. Y. Sonin.
The interval of definition is \((0, \infty )\), the weight function \(w(x)=e^{-x}x^a\) and, for the orthonormal case, the leading coefficients are given by \(k_n=(-1)^n/n!\), \(\tilde {h}_n = \Gamma (n+a+1)/n!\) and \(r_n=-n(n+a)\). The parameter \(a\) should satisfy \(a > -1\).
The function call HermiteP(n,x) will return the \(n\)th Hermite polynomial if \(n\) is a non-negative integer; otherwise the result will involve the original operator HermiteP or on graphical interfaces \(H_n(x)\) will be output.
The interval of definition is \((-\infty , +\infty )\), the weight function \(w(x)=e^{-x^2}\) and, for the orthonormal case, the leading coefficients are given by \(k_n=2^n\), \(\tilde {h}_n = \sqrt {\pi }2^nn!\) and \(r_n=0\).
FibonacciP(n,x) returns the \(n\)th Fibonacci polynomial in the variable \(x\). If \(n\) is an integer, it will be evaluated using the recursive definition: \[F_0(x) = 0;\qquad F_1(x) = 1; \qquad F_n(x) = x F_{n-1}(x) + F_{n-2}(x)\,.\]
The recursion is, of course, optimised as a simple loop to avoid repeated computation of lower-order polynomials.
Euler numbers are computed by the unary operator Euler; the call Euler(n) returns the \(n\)th Euler number; all the odd Euler numbers are zero. The computation is derived directly from Pascal’s triangle of binomial coefficients.
The Euler numbers and polynomials have the following generating functions:
\[\frac {2e^t} {1+e^{2t}} = \sum _{n=0}^{\infty }\frac {E_nt^n}{n!},\qquad \quad \frac {e^{x t}} {1+e^t} = \sum _{n=0}^{\infty }\frac {E_n(x)t^n}{n!}\] respectively. Thus \(E_0=1\) and \(E_1=0\). Furthermore the numbers and polynomials are related by the equations: \[ E_n = 2^nE_n\left (\frac {1}{2}\right ),\qquad \quad E_n(x) = \sum _{k=0}^n \binom {n}{k} \frac {E_k}{2^k} \left (x-\frac {1}{2}\right )^{n-k}.\]
The Euler polynomials are evaluated for non-negative integer \(n\) by using the summation immediately above.
The call Bernoulli(n) evaluates to the \(n\)th Bernoulli number; all of the odd Bernoulli numbers, except Bernoulli(1), are zero.
The algorithms for Bernoulli numbers used are based upon those by Herbert Wilf, presented by Sandra Fillebrown [Fil92]. If the rounded switch is off, the algorithms are exactly those; if it is on, some further rounding may be done to prevent computation of redundant digits. Hence, these functions are particularly fast when used to approximate the Bernoulli numbers in rounded mode.
The Bernoulli numbers and polynomials have the following generating functions: \[ \frac {t} {e^t-1} = \sum _{n=0}^{\infty }\frac {B_nt^n}{n!},\qquad \quad \frac {te^{x t}} {e^t -1} = \sum _{n=0}^{\infty }\frac {B_n(x)t^n}{n!}\] respectively. Thus \(B_0=1\) and \(B_1=-\frac {1}{2}\). Furthermore the numbers and polynomials are related by the equations: \[ B_n= B_n(0),\qquad \quad B_n(x) = \sum _{k=0}^n \binom {n}{k} B_kx^{n-k}.\]
The Bernoulli polynomials are evaluated for non-negative integer \(n\) by using the summation immediately above.
Both the Bernoulli and Euler numbers and polynomials may also be calculated directly by expanding the corresponding generating function as a power series in \(t\) using either the TPS or TAYLOR package, extracting the \(n\)th term and multiplying by \(n!\). The use of the TPS package is probably preferable here as the series for the generating function is extendible and need only be calculated once; it will be extended automatically if higher order numbers or polynomials are required.
The following procedures compute sets of functions e.g. to be used for approximation. All procedures have two parameters, the expression to be used as \(variable\) (an identifier in most cases) and the order of the desired system. The functions are not scaled to a specific interval, but the \(variable\) can be accompanied by a scale factor and/or a translation in order to map the generic interval of orthogonality to another (e.g. \((x- 1/2 ) * 2 pi\)). The result is a function list with ascending order, such that the first element is the function of order zero and (for the polynomial systems) the function of order \(n\) is the \(n+1\)-th element.
monomial_base(x,n) {1,x,...,x**n}
trigonometric_base(x,n) {1,sin x,cos x,
sin(2x),cos(2x)...}
Bernstein_base(x,n) Bernstein polynomials
Legendre_base(x,n,a,b) Legendre polynomials
Laguerre_base(x,n,a) Laguerre polynomials
Hermite_base(x,n) Hermite polynomials
Chebyshev_base_T(x,n) Chebyshev polynomials
of the first kind
Chebyshev_base_U(x,n) Chebyshev polynomials
of the second kind
Gegenbauer_base(x,n,a) Gegenbauer polynomials
Example:
Bernstein_base(x,5);
5 4 3 2
{ - X + 5*X - 10*X + 10*X - 5*X + 1,
4 3 2
5*X*(X - 4*X + 6*X - 4*X + 1),
2 3 2
10*X *( - X + 3*X - 3*X + 1),
3 2
10*X *(X - 2*X + 1),
4
5*X *( - X + 1),
5
X }
The contributions of Kerry Gaskell, Matthew Rebbeck, Lisa Temme, Stephen Scowcroft and David Hobbs (all students from the University of Bath on placement in ZIB Berlin for one year) were very helpful to augment the package. The advice of René Grognard (CSIRO, Australia) for the development of the module for Clebsch-Gordan and 3j, 6j symbols and the module for spherical and solid harmonics was very much appreciated.
Special Functions
| Function | Operator |
| \(\mathrm {Si}(z)\) | Si(z) |
| \(\mathrm {Si}(z)-\pi /2\) | s_i(z) |
| \(\mathrm {Ci}(z)\) | Ci(z) |
| \(\mathrm {Shi}(z)\) | Shi(z) |
| \(\mathrm {Chi}(z)\) | Chi(z) |
| \(\mathrm {erf}(z)\) | Erf(z) |
| \(1-\mathrm {erf}(z)\) | erfc(z) |
| \(\mathrm {Ei}(z)\) | Ei(z) |
| \(\mathrm {Ei}(\log (z))\) | Li(z) |
| \(C(x)\) | Fresnel_C(x) |
| \(S(x)\) | Fresnel_S(x) |
| \(B(a,b)\) | Beta(a,b) |
| \(\Gamma (a)\) | Gamma(a) |
| normalized incomplete Beta \(I_{x}(a,b)=\frac {\textstyle B_{x}(a,b)}{\textstyle B(a,b)}\) | iBeta(a,b,x) |
| normalized incomplete Gamma \(P(a,z)=\frac {\textstyle \gamma (a,z)}{\textstyle \Gamma (a)}\) | iGamma(a,z) |
| incomplete Gamma \(\gamma (a,z)\) | m_gamma(a,z) |
| \((a)_k\) | Pochhammer(a,k) |
| \(\psi (z)\) | Psi(z) |
| \(\psi ^{(n)}(z)\) | Polygamma(n,z) |
| \(J_\nu (z)\) | BesselJ(nu,z) |
| \(Y_\nu (z)\) | BesselY(nu,z) |
| \(I_\nu (z)\) | BesselI(nu,z) |
| \(K_\nu (z)\) | BesselK(nu,z) |
| \(H^{(1)}_n(z)\) | Hankel1(n,z) |
| \(H^{(2)}_n(z)\) | Hankel2(n,z) |
| Function | Operator |
| \(\mathrm {Ai}(z)\) | Airy_Ai(z) |
| \(\mathrm {Bi}(z)\) | Airy_Bi(z) |
| \(\mathrm {Ai}'(z)\) | Airy_Aiprime(z) |
| \(\mathrm {Bi}'(z)\) | Airy_Biprime(z) |
| \(\mathbf {H}_{\nu }(z)\) | StruveH(nu,z) |
| \(\mathbf {L}_{\nu }(z)\) | StruveL(nu,z) |
| \(s_{a,b}(z)\) | Lommel1(a,b,z) |
| \(S_{a,b}(z)\) | Lommel2(a,b,z) |
| \(M(a, b, z)\) or \(_1F_1(a, b; z)\) or \(\Phi (a, b; z)\) | KummerM(a,b,z) |
| \(U(a, b, z)\) or \(z^{-a}{_2F_0(a, b; z)}\) or \(\Psi (a, b; z)\) | KummerU(a,b,z) |
| Expression in Kummer_M | WhittakerM(kappa,mu,z) |
| Expression in Kummer_U | WhittakerW(kappa,mu,z) |
| Riemann’s \(\zeta (z)\) | zeta(z) |
| Lambert \(\omega (z)\) | Lambert_W(z) |
| \(\mathrm {Li}_2(z)\) | dilog(z) |
| \(\mathrm {Li}_{n}(z)\) | Polylog(n,z) |
| Lerch’s transcendent \(\Phi (z,s,a)\) | Lerch_Phi(z,s,a) |
\( \newcommand {\ThreeJSymbol }{\begin {pmatrix} j_{1} & j_{2} & j_{3} \\ m_{1} & m_{2} & m_{3} \end {pmatrix}} \newcommand {\SixJSymbol }{\begin {Bmatrix} j_{1} & j_{2} & j_{3} \\ m_{1} & m_{2} & m_{3} \end {Bmatrix}} \)
| Function | Operator |
| \(Y_n^{m}(x,y,z,r2)\) | SolidHarmonicY(n,m,x,y,z,r2) |
| \(Y_n^{m}(\theta ,\phi )\) | SphericalHarmonicY(n,m,theta,phi) |
| \(\displaystyle \ThreeJSymbol \) | ThreeJSymbol({j1,m1},{j2,m2}, |
| \(( j_1m_1j_2m_2 \mid j_1j_2j_3 - m_3 )\) | Clebsch_Gordan({j1,m1},{j2,m2}, |
| \(\displaystyle \SixJSymbol \) | SixJSymbol({j1,j2,j3},{l1,l2,l3}) |
| Function | Operator |
| Fibonacci Polynomials \(F_{n}(x)\) | FibonacciP(n,x) |
| \(B_n(x)\) | BernoulliP(n,x) |
| \(E_n(x)\) | EulerP(n,x) |
| \(H_n(x)\) | HermiteP(n,x) |
| \(L_n(x)\) | LaguerreP(n,x) |
| Generalised Laguerre \(L_n^{(m)}(x)\) | LaguerreP(n,m,x) |
| \(P_n(x)\) | LegendreP(n,x) |
| Associated Legendre \(P_n^{(m)}(x)\) | LegendreP(n,m,x) |
| \(U_n(x)\) | ChebyshevU(n,x) |
| \(T_n(x)\) | ChebyshevT(n,x) |
| \(C_n^{(\alpha )}(x)\) | GegenbauerP(n,alpha,x) |
| \(P_n^{(\alpha ,\beta )} (x)\) | JacobiP(n,alpha,beta,x) |
Well-known Numbers and Reserved Constants
| Function | Operator |
| \(\displaystyle \binom {n}{m}\) | Binomial(n,m) |
| Fibonacci Numbers \(F_{n}\) | Fibonacci(n) |
| \(\mathrm {s}_n^m\) | Stirling1(n,m) |
| \(\mathrm {S}_n^m\) | Stirling2(n,m) |
| Bernoulli(\(n\)) or \( B_n \) | Bernoulli(n) |
| Euler(\(n\)) or \( E_n \) | Euler(n) |
| Motzkin(\(n\)) or \( M_n \) | Motzkin(\(n\)) |
| Constant | REDUCE name |
| Square root of \((-1)\) | i |
| \(\pi \) | pi |
| Base of natural logarithms | e |
| Euler’s \(\gamma \) constant | Euler_gamma |
| Catalan’s constant | Catalan |
| Khinchin’s constant | Khinchin |
| Golden ratio | Golden_ratio |
| Up | Next | Prev | PrevTail | Front |
Hosted by