SPECFN: Package for special functions

REDUCE

16.68 SPECFN: Package for special functions

This special function package is separated into two portions to make it easier to handle. The packages are called SPECFN and SPECFN2. The ﬁrst one is more general in nature, whereas the second is devoted to special special functions. Additional documentation and examples can be found in the ﬁles specfn.tex, specfn.tst and specfn2.tst in the packages/specfn directory.

Author: Chris Cannam, with contributions from Winfried Neun, Herbert Melenk, Victor Adamchik, Francis Wright, Alan Barnes and several others.

16.68.1 Special Functions: Introduction

The package SPECFN is designed to provide algebraic and numeric manipulations of many common special functions, namely:

The Exponential Integral, Sine & Cosine Integrals;
The Hyperbolic Sine & Cosine Integrals;
The Fresnel Integrals & Error function;
The Gamma function;
The Beta function;
The psi function & its derivatives;
The Bessel functions J and Y of the ﬁrst and second kinds;
The modiﬁed Bessel functions I and K;
The Hankel functions H⁽¹⁾ and H⁽²⁾;
The Airy functions;
The Kummer hypergeometric functions M and U;
The Struve, Lommel and Whittaker functions;
The Riemann Zeta function;
The Dilog function;
The Polylog and Lerch Phi functions;
Lambert’s W function;
Associated Legendre Functions (Spherical and Solid Harmonics);
3j and 6j symbols, Clebsch-Gordan coeﬃcients;
Jacobi’s Elliptic Functions;
Elliptic Integrals;
Jacobi’s Theta Functions;
Nome of Elliptic Functions;
Stirling Numbers;
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].

There is a nice collection of plot calls for special functions in the ﬁle specplot.tst in the subfolder plot of the packages folder. These examples will reproduce a number of well-known pictures from [AS72].

16.68.2 Polynomial Functions: Introduction

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;

16.68.2.1 Orthogonal Polynomial Functions

The polynomial function sets available are:

Hermite Polynomials;
Legendre Polynomials;
Laguerre Polynomials;
Chebyshev Polynomials;
Jacobi Polynomials;
Gegenbauer Polynomials;

16.68.3 Simpliﬁcation and Approximation

All of the operators supported by this package have certain algebraic simpliﬁcation 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 simpliﬁcations within the special function package are deﬁned in the form of a REDUCE ruleset. Therefore, in order to get a quick insight into the simpliﬁcation 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 diﬀerent (hopefully better) results for the algebraic simpliﬁcations 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 oﬀ in order to inhibit the storage of approximated values. The switch is on by default.

16.68.4 Integral Functions

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 deﬁned by:

∫ z erf(z) = √2-- e− t2dt π 0 2 ∫ ∞ 2 erfc(z) = √--- e−t dt = 1 − erf(z ) ∫ π z z (π- 2) C (z) = cos 2 t dt ∫0z ( ) S(z) = sin π-t2 dt 0 2

respectively.

The exponential and related integrals are deﬁned by the following:

∫ ∞ e− t Ei(z) = e− z -----dt ∫ z t + z z-dt-- Li(z) = logt ∫ 0z Si(z) = sintdt 0 t ∫ ∞ sint π si(z) = − ----dt = Si(z) − -- ∫ z∞ t ∫ z 2 Ci(z) = − costdt = cos-t −-1-dt + log z + γ z t 0 t ∫ zsinht Shi(z) = -----dt ∫ 0 t zcosh-t −-1 Chi (z) = t dt + log z + γ 0

where γ is Euler’s constant (Euler_gamma).

The deﬁnitions 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 suﬀer 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.

Li(z) is simpliﬁed to Ei(ln(z)) .

16.68.5 The Γ Function and Related Functions

16.68.5.1 The Γ Function

This is represented by the unary operator Gamma. The Gamma function is deﬁned by the integral:

∫ ∞ − ta−1 Γ (a) = e t dt. 0

Initial transformations applied with ROUNDED oﬀ are: Γ(n) for integral n is computed, Γ(n + 1∕2) for integral n is rewritten to an expression in √ --
π

, Γ(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 Γ(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(Γ), with a scaling factor obtained from the Pochhammer symbols.

An expression for Γ′(z) in terms of Γ and ψ is included.

16.68.5.2 Incomplete Gamma Functions

The (unnormalised) incomplete gamma function is provided by the binary function m_gamma. In the literature it is normally represented as γ(a,z) and is deﬁned by

∫ z γ(a,z) = e−tta− 1dt. 0

The normalised incomplete gamma function P(a,z) is provided by the binary function igamma and is deﬁned as

P (a,z) = γ(a,z). Γ (a)

16.68.5.3 The Beta Functions

The binary function B(a,b) is related to the Γ function[AS72] and is deﬁned by

∫ 1 Γ (a)Γ (b) B (a, b) = ta(1 − t)bdt = ---------. 0 Γ (a + b)

It is represented by the binary function Beta.

The unnormalised and nomalised incomplete Beta funtions are deﬁned by

∫ x a b Bx(a,b) = 0 t (1 − t)dt, Ix(a,b) = Bx-(a,b) B (a,b)

respectively. The normalised one is represented by the ternary function ibeta(a,b,x).

16.68.5.4 The Digamma Function, ψ

This is represented by the unary operator psi. It is deﬁned as the logarithmic derivative of the Γ function:

Γ-′(z) ψ (z) = Γ (z).

Initial transformations for ψ are applied on a similar basis to those for Γ; where possible, ψ(x) is rewritten in terms of ψ(1) and ψ(), and expressions with negative arguments are rewritten to have positive ones.

The algorithm for numerical evaluation of ψ is based upon an asymptotic series, with a suitable scaling.

Relations for the derivative and integral of ψ are included.

16.68.5.5 The Polygamma Functions, ψ⁽ⁿ⁾

The nth derivative of the ψ function is represented by the binary operator Polygamma, whose ﬁrst argument is n.

Initial manipulations on ψ⁽ⁿ⁾ are few; where the second argument is 1 or 3∕2, the expression is rewritten to one involving the Riemann ζ function, and when the ﬁrst is zero it is rewritten to ψ; 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 Reﬂection Formula is used, introducing a term in the nth derivative of cot(zπ).

Simple relations for derivatives and integrals are provided.

16.68.6 Bessel Functions

Support is provided for the Bessel functions J and Y , the modiﬁed Bessel functions I and K, and the Hankel functions of the ﬁrst and second kinds. The relevant operators are, respectively, BesselJ, BesselY, BesselI, BesselK, Hankel1 and Hankel2, which are all binary operators.

The Bessel functions J_ν(z) and Y _ν(z) are solutions of the Bessel equation:

2 z2 d-w-+ z dw-+ (z2 − ν2 )w = 0. dz2 dz

Bessel’s function of the ﬁrst kind, J_ν(z), has the series expansion:

( ) ν∑∞ 2k J (z) = z- (− 1)k----(z∕2)-----. ν 2 k!Γ (ν + k + 1) k=0

Bessel’s function of the second kind, Y _ν(z), (for non-integral ν) is deﬁned by:

Yν(z) = Jν(z)-cos(νπ)-−-J−ν(z) sin(ν π)

or by its limiting value:

|| n || Y ν(z ) = 1-∂Jν(z)| + (−-1)-∂J-ν(z-)| . π ∂ν |ν=n π ∂ν |ν= −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 →∞:

∘ --- (1) 2 ( ( νπ π)) Hν (z) ∼ ---exp i z − ---− -- , ∘ πz- 2 4 (2) 2-- ( ( ν-π π-)) Hν (z) ∼ πz exp − i z − 2 − 4 .

The modiﬁed Bessel functions I_ν(z) and K_ν(z) are solutions of the modiﬁed Bessel equation:

d2w dw z2 ----+ z ---− (z2 + ν2 )w = 0. dz2 dz

Since they may be obtained by replacing z by ±iz the modiﬁed Bessel functions are sometimes called Bessel functions of imaginary argument. I_ν(z) has the series expansion:

( )ν∑∞ 2k Iν(z) = z- ----(z∕2)------, 2 k!Γ (ν + k + 1 ) k=0

whereas K_ν(z) is distinguished by its asymptotic behaviour:

∘ --- π − z Kν(z) ∼ 2ze

as z →∞. For more information, see the DLMF:NIST chapters on Hankel & Bessel functions and Modiﬁed 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 ﬁrst argument are transformed to have positive ﬁrst 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 ﬁrst argument are carried out;
derivatives of Hankel functions are carried out.

Also, if the COMPLEX switch is on and ROUNDED is oﬀ, 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 ﬁrst, 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.

16.68.7 Airy 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 diﬀerential equation:

2 d-w- = zw. dz2

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:

1 √ -[ ] 2 2 Ai (z) = -- z I−1∕3(ζ) − I1∕3(ζ) , where ζ = --z3. 3 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.

16.68.8 Hypergeometric and Other Functions

This package also provides some support for other functions, in the form of algebraic simpliﬁcations:

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 diﬀerentiate with respect to the second argument. These functions with arguments ν and x are solutions of the diﬀerential equation: $2 ( 2) ν−1 d-w- + 1-dw- + 1 − ν-- w = √--(z∕2)------. dx2 x dx x2 πΓ (ν + 1∕2)$
The Lommel functions of the ﬁrst and second kinds, through the ternary operators Lommel1 and Lommel2 with arguments ν, μ and x may be considered generalisations of the Struve functions satisfying the diﬀerential equation: $d2w 1 dw ( ν2) ---2 + ----- + 1 − -2- w = zμ−1. dx x dx x$ Manipulations are provided to handle special cases and simplify where appropriate.
The Kummer conﬂuent hypergeometric functions M and U (the hypergeometric ₁F₁ or Φ, and z^−a₂F₀ or Ψ, respectively), represented by the ternary operators KummerM and KummerU with arguments a, b and x, are solutions of the diﬀerential equation: $d2w dw --2-+ (b − x) ---− aw = 0. dx dx$ There are manipulations for special cases and simpliﬁcations, 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 κ, μ and x. They satisfy the Whittaker diﬀerential equation: $2 ( 2 ) d-W--+ 1-−-4μ--+ κ-− 1- W = 0, dx2 4x2 x 4$ which is obtained from the Kummer diﬀerential equation via the substituions $W = ez∕2zμ+1 ∕2w, κ = b∕2 − a μ = (b − 1 )∕2.$ The Whittaker M and W functions with non-numeric arguments are simpliﬁed to expressions involving the Kummer M and U functions respectively.

16.68.9 The Riemann Zeta Function

This is represented by the unary operator Zeta and deﬁned by the formula:

∑∞ 1 ζ(s) = --. n=1 ns

With ROUNDED oﬀ, ζ(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 ζ is only carried out if the argument is real. The algorithms used for ζ are: for odd integral arguments, an expression relating ζ(n) with ψⁿ⁻¹(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 ﬁrst few primes in the standard over-all-primes expansion, and then continue with the deﬁning 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 diﬀerentiation or integration of ζ.

16.68.10 Polylogarithm and Related Functions

The dilogarithm function Li₂(z) is deﬁned by

∫ ∑∞ zn zlog(1 − t) Li2(z) ≡ n2-= − -----t----dt n=1 0

and represented by the unary function dilog.

The polylogarithm function Li_s(z) is deﬁned by

∑∞ n ∫ ∞ s− 1 Lis(z) ≡ z--= --z-- -t----dt. n=1 ns Γ (s) 0 et − z

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 ζ(s). For s = 1, the polylogarithm reduces to the elementary function: − log(1 − t).

Lerch’s transcendent or Lerch Phi function is deﬁned by

∑∞ zn Φ (z, s,a) = --------. n=0 (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: zLi_s(z) = Φ(z,s, 1).

16.68.11 Lambert’s W Function

Lambert’s function ω(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 ω(x) is only real-valued in the interval (−1∕e,∞). In the interval (−1∕e, 0), it is double-valued with a branch point at the point (-1/e, -1) where ω′(x) is singular. The positive branch is deﬁned on the interval (−1∕e,∞) where it is monotonically increasing with ω(x) > −1. The negative branch is deﬁned on the interval (−1∕e, 0) where it is monotonically decreasing with ω(x) < −1.

Simpliﬁcation rules for ω(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 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.

16.68.12 Spherical and Solid Harmonics

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 deﬁnition:

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 ≥ 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≠0.

S(n,m,x,y,z) is proportional to

r^n * Legendre(n,m,cos theta) * exp(i*phi)

where r² = x² + y² + z².

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 ≥ 0,−n ≤ m ≤ n form a complete orthogonal basis on it, i.e. ⟨n,m|n′,m′⟩ = δ_n,n′δ_m,m′ using ⟨…|…⟩ to designate the scalar product of functions over the spherical surface.

The coeﬃcients 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 deﬁnition.

It is also clear that all diﬀerentiations 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),…, (n,m).

Numerically these recursions are stable.

For m < 0 one has:-

S (n,m, x,y,z ) = (− 1)mS (n,− m, x,− y,z).

16.68.13 3j symbols and Clebsch-Gordan Coeﬃcients

The operators ThreeJSymbol and Clebsch_Gordan are deﬁned 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});

16.68.14 6j symbols

The operator SixJSymbol is deﬁned as in [LB68] or [Edm57] and expects two lists of values {j₁,j₂,j₃} and {l₁,l₂,l₃} 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.

16.68.15 Elliptic Functions

The implementation of the functions in this and the next two subsections have, in 2019, been substantially revised by Alan Barnes. This is to bring the notation more into line with standard (British) texts such as Whittaker & Watson [WW69] and Lawden [Law89] and also to correct a number of errors and omissions. These changes and additions will be itemised in the relevant sections below.

The following functions have been implemented:

The 12 Jacobi Functions
Jacobi’s Amplitude Function
Arithmetic Geometric Mean
Descending Landen Transformation

16.68.15.1 Jacobi Elliptic Functions

The following Jacobi functions are available:-

jacobisn(u,k)
jacobidn(u,k)
jacobicn(u,k)
jacobicd(u,k)
jacobisd(u,k)
jacobind(u,k)
jacobidc(u,k)
jacobinc(u,k)
jacobisc(u,k)
jacobins(u,k)
jacobids(u,k)
jacobics(u,k)

These diﬀer somewhat from the originals implemented by Lisa Temme in that the second argument is now the modulus (usually denoted by k in most texts rather than its square m). The notation for the most part follows Lawden [Law89]. The last nine Jacobi functions are related to the three basic ones: jacobisn(u,k), jacobicn(u,k) and jacobidn(u,k) and use Glaisher’s notation. For example

ns(x,k) = ---1----, cs(x, k) = cn(u,k-), cd(x, k) = cn(u,k-). sn(u, k) sn(u,k) dn(u, k)

Extended rule lists are provided for diﬀerentiation of these functions with respect to either argument, to implement the standard addition formulae, argument shifts by multiples of the two quarter-periods K and iK′ and ﬁnally Jacobi’s transformation for a purely imaginary ﬁrst argument.

When their arguments are purely numerical, these functions will be evaluated numerically if the rounded switch is used. For complex arguments it is also better if the complex switch is on.

16.68.15.2 Jacobi Amplitude Function

The amplitude of u can be evaluated using the jacobiam(u,k) command. A rule list is provided for diﬀerentiation of this functions with respect to either argument.

16.68.15.3 Arithmetic Geometric Mean (AGM)

A procedure to evaluate the AGM of initial values $a0,b0,c0$ exists as
AGM_function( $a0,b0,c0$ ) and will return
{N,AGM,{a_N,…,a₀},{b_N,…,b₀},{c_N,…,c₀}}, where N is the number of steps to compute the AGM to the desired accuracy.

To determine the Elliptic Integrals K(m), E(m) we use initial values $a0 = 1$ ; $√ ------- b0 = 1 − k2$ ; $c0 = k$ .

This procedure and the following one are primarily intended for use in the numerical evaluation of the various elliptic functions and integrals rather than for direct use by users.

16.68.15.4 Descending Landen Transformation

The procedure to evaluate the Descending Landen Transformation of ϕ and α uses the following equations:

(1 + sin α_n+1)(1 + cos α_n) = 2	where α_n+1 < α_n,
tan(ϕ_n+1 − ϕ_n) = cos α_n tan ϕ_n	where ϕ_n+1 > ϕ_n.

It can be called using landentrans(ϕ₀, α₀) and will return
{{ϕ₀,…,ϕ_n},{α₀,…,α_n}}.

16.68.16 Elliptic Integrals

The following functions have been implemented:

Complete & Incomplete Elliptic Integrals of the First Kind
Complete & Incomplete Elliptic Integrals of the Second Kind
Jacobi’s Zeta Function

These again diﬀer somewhat from the originals implemented by Lisa Temme as the second argument is now the modulus k rather that its square. Also in the original implementation there was some confusion between Legendre’s form and Jacobi’s form of the incomplete elliptic integrals of the second kind; E(u,k) denoted the ﬁrst in numerical evaluations and the second in the derivative formulae for the Jacobi elliptic functions with respect to their second argument. This confusion was perhaps understandable as in the literature some authors use the notation E(u,k) for the Legendre form and others for Jacobi’s form.

To bring the notation more into line with that in the NIST Digital Library of Mathematical Functions and avoid any possible confusion, E(u,k) is used for the Legendre form and ℰ(u,k) for Jacobi’s form. This diﬀers from the 2019 version of this section which followed Lawden [Law89], where the notation D(ϕ,k) and E(u,k) were used for the Legendre and Jacobi forms respectively.

A number of rule lists have been provided to implement, where appropriate, derivatives of these functions, addition rules and periodicity and quasi-periodicity properties and to provide simpliﬁcations for special values of the arguments.

16.68.16.1 Elliptic F

The Elliptic F function can be used as EllipticF(phi,k) and will return the value of the Incomplete Elliptic Integral of the First Kind:

∫ ϕ F(ϕ,k) = (1 − k2 sin2 𝜃)−1∕2d𝜃. 0

16.68.16.2 Elliptic K

The Elliptic K function can be used as EllipticK(k) and will return the value of the Complete Elliptic Integral of the First Kind:

∫ π∕2 2 2 −1∕2 K(k ) = F (π∕2,k ) = (1 − k sin 𝜃) d𝜃. 0

This is one of the quarter periods of the Jacobi elliptic functions and is often used in the calculation of other elliptic functions.

The complementary Elliptic K′ function can be used as EllipticK!′(k) and will return the value

′ √ ------- K (k ) = K( 1 − k2)

which is the other quarter-period of the Jacobi elliptic functions.

16.68.16.3 Elliptic E

The Elliptic E function comes with either one or two arguments; used with two arguments as EllipticE(u,k) it will return the value of Legendre’s form of the Incomplete Elliptic Integral of the Second Kind:

∫ ϕ∘ -----2---2-- E(ϕ, k) = 1 − k sin 𝜃d 𝜃. 0

When called with one argument EllipticE(k) will return the value of the Complete Elliptic Integral of the Second Kind:

∫ π∕2∘ ------------ E (k) = E(π∕2, k) = 1 − k2 sin2 𝜃d𝜃. 0

The complementary Elliptic E′ function can be used as EllipticE!′(k) and will return the value

√ ------- E (k′) = E ( 1 − k2 ).

16.68.16.4 Jacobi E

The Jacobi E function can be used as JacobiE(u,k); it will return the value of Jacobi’s form of the Incomplete Elliptic Integral of the Second Kind:

∫ u ℰ(u,k) = dn2(v,k )dv. 0

The relationship between the two forms of incomplete elliptic integrals can be expressed as

ℰ (u, k) = E(am (u),k ).

Note that

∫ K(k) 2 E(k ) = ℰ (K(k ),k ) = dn (v,k)dv. 0

Display of the Jacobi E function in REDUCE causes some minor diﬃculty. If the GUI supports calligraphic characters there is no problem and it is rendered as ℰ(u,k). However, this is currently not the case with the CSL GUI and so it is rendered as 𝜖. The uppercase version cannot reasonably be used as it is indistinguishable from uppercase E which is used for Elliptic E. In non-GUI interfaces the Jacobi E function is rendered as Epsilon.

16.68.16.5 Jacobi’s Zeta Function

This can be called as JacobiZeta(u,k) and refers to Jacobi’s (elliptic) Zeta function Z(u,k) whereas the operator Zeta will invoke Riemann’s ζ function.

16.68.17 Jacobi Theta Functions

These theta functions diﬀer from those originally deﬁned by Lisa Temme in a number of respects. Firstly four separate functions of two arguments are deﬁned:

elliptictheta1(u,q) 𝜗₁(u,q)
elliptictheta2(u,q) 𝜗₂(u,q)
elliptictheta3(u,q) 𝜗₃(u,q)
ellipticthetas(u,q) 𝜗₄(u,q)

rather than a single function with three arguments (with the ﬁrst argument taking integer values in the range 1 to 4). Secondly the periods are 2π, 2π,π and π respectively (NOT 4K, 4K, 2K and 2K). Thirdly the second argument is the nome q where |q| < 1 and hence the quasi-period is −i log q.

elliptictheta1 and elliptictheta2 are consequently multi-valued owing to the appearance of q^1∕4 in their deﬁning expansions. elliptictheta3 and elliptictheta4 are, however, single-valued functions of q.

Regarded as functions of the ‘parameter’ (usually denoted by τ where q = exp(iπτ)), elliptictheta1 and elliptictheta2 are single-valued functions of τ.

Recall that in order that |q| < 1, the imaginary part of τ must be positive. Note also in this case q^1∕4 is interpreted as exp(iπτ∕4) rather than the principal value of q^1∕4. Thus τ, 2 + τ, 4 + τ and 6 + τ produce four diﬀerent values of both elliptictheta1 and elliptictheta2 although they all produce the same nome q.

Utilising the periodicity and quasi-periodicity of the theta functions some generalised shift rules are implemented to shift their ﬁrst argument into the base period parallelogram with vertices

(π∕2,π τ∕2), (− π∕2, πτ∕2 ), (− π ∕2,− πτ∕2), (π∕2,− πτ∕2 ).

Together with the relation 𝜗₁(0,q) = 0, these shift rules serve to simplify all four theta functions to zero when appropriate.

When the switch rounded is on and the arguments are purely numerical, the theta functions are evaluated numerically. For complex numerical arguments the switch complexmust also be on.

The numerical code always returns the principal value of elliptictheta1 and elliptictheta2 (that is the one obtained by taking the principal value of q^1∕4). Thus the negative real-axis of q is a branch cut with the cut line regarded as belonging to the second rather than the third quadrant.

The series for the theta functions are fairly rapidly convergent due to the quadratic growth of the exponents of q – except for values of q for which |q| is near to 1. In such cases the direct algorithm would suﬀer from slow convergence and rounding errors. For such values of |q|, Jacobi’s transformation τ′ = −1∕τ can be used to produce a smaller value of the nome and so increase the rate of convergence. When ℜq < 0, the Jacobi transformation is preceded by the modular transformation τ′ = 1 + τ, which has the eﬀect of multiplying q by −1, so that the new nome has a non-negative real part.

This works very well for real values of q, but for complex values the gains are somewhat smaller. For a nome q whose modulus is near to 1, the Jacobi transformation produces a nome q′ for which |q′|≈|q|^α² where α = π∕| arg(q)|≥ 2. Thus the worst case occurs for values of the nome q near to ±i where α ≈ 2.

16.68.17.1 Some Numerical Utility Functions

Five utility functions are provided:

nome2mod(q)
nome2mod!′(q)
nome2!K(q)
nome2!K!′(q)
nome(k)

These are only operative when the switch rounded is on and their argument is numerical. The ﬁrst pair relate the nome q of the theta functions with the moduli k and k′ = √1-−--k2- of the associated Jacobi elliptic functions.

The second pair return the quarter-periods K and K′ respectively of the Jacobi elliptic functions associated with the nome q.

Finally, nome(k) returns the nome q associated with the modulus k of a Jacobi elliptic function and is essentially the inverse of nome2mod.

16.68.18 Stirling Numbers

The Stirling numbers of the ﬁrst and second kind are computed by calling the binary operators Stirling1 and Stirling2 respectively.

Stirling numbers of the ﬁrst kind have the generating function:

∑n smn xm = (x − n + 1)n m=0

where (x−n + 1)_n is the Pochhammer symbol. This provides a convenient way of calculating these Stirling numbers by extracting coeﬃcients of the polynomial obtained by evaluating the Pochhammer symbol. REDUCE however uses an explicit summation.

Stirling numbers of the second kind are deﬁned by the formula:

1 ∑m (n ) Smn = --- (− 1)m −k kn. m! k=0 k

REDUCE uses this explicit summation to evaluate Stirling numbers of the second kind.

16.68.19 Constants

The following well-known constants are deﬁned 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 −ψ(1);
Catalan : Catalan’s constant;
Khinchin : Khinchin’s constant, deﬁned in [Khi64] (which takes a lot of time to compute);
Golden_Ratio :

16.68.20 Orthogonal Polynomials

All the polynomials in this section take two or more parameters; the ﬁrst 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 deﬁnitions appropriate to all the sets of orthogonal polynomials in the following subsections are listed.

A set of polynomials {p_n(x)},n = 0, 1,… are said to be orthogonal on open interval (a,b) (where a and/or b may be inﬁnite) with positive weight function w(x) if

∫ b pn(x )pm (x)w (x)dx = 0 when m ⁄= n. a

This deﬁnes each polynomial p_n(x) up to a constant factor c_n which is usually ﬁxed by normalisation. If these factors are chosen so that

∫ b 2 ∘ --- hn = a (pn(x)) w(x )dx = 1 i.e. cn = hn

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≠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 coeﬃcient 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 ‘ﬁrst moment’ h_n of the nth polynomial deﬁned by:

∫ ˜ b 2 hn = x (pn(x))w (x)dx a

and the ratio

˜k rn = -n- where pn(x) = knxn + ˜knxn −1... kn

These quantities may be used in recurrence relations when generating the polynomials.

16.68.20.1 Legendre Polynomials

The function call LegendreP(n,x) will return the nth 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 deﬁnition is (−1, 1), the weight function w(x) = 1 and, for the orthonormal case, the leading coeﬃcients are given by k_n = 2ⁿ( 1
2 )_n∕n! where ( 1
2 )_n is the Pochhammer symbol. Also h_n = --2-
2n+1 and r_n = 0.

16.68.20.2 Associated Legendre Functions

The function call LegendreP(n,m,x) will return the nth associated Legendre function if n and m are integers with 0 ≤ m ≤ n; otherwise the result will involve the original operator LegendreP or on graphical interfaces P_n^(m)(x) will be output.

They are deﬁned by

dmPn (x) P(nm)(x) = (− 1)m(1 − x2)m ∕2-----m---; dx

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 ﬁxed 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:

∫ 1 ( ) P(nm)(x) 2 dx = hn = ---2(l +-m-)!----. −1 (2l + 1)!(l − m )!

16.68.20.3 Chebyshev Polynomials

The function call ChebyshevT(n,x) will return the nth Chebyshev polynomial of the ﬁrst 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 deﬁnition is (−1, 1), the weight function w(x) = (1 − x²)^−1∕2 and, for the orthonormal case, the leading coeﬃcients are given by k_n = 2ⁿ⁻¹ for n > 0; k₀ = 1. Also h_n = π∕2 for n > 0; h₀ = π and r_n = 0.

The function call ChebyshevU(n,x) will return the nth 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 deﬁnition is (−1, 1), the weight function w(x) = (1 − x²)^−1∕2 and, for the orthonormal case, the leading coeﬃcients are given by k_n = 2ⁿ, h_n = π∕2 and r_n = 0.

16.68.20.4 Gegenbauer Polynomials

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 deﬁnition is (−1, 1), the weight function w(x) = (1 −x²)^a−1∕2 and, for the orthonormal case, the leading coeﬃcients are given by k_n = 2ⁿ(a)_n∕n! where (a)_n is the Pochhammer symbol. The parameter a should satisfy a > −1∕2, a≠0. Also

˜ 21−-2aπ-Γ (n-+-2a) hn = (n + a)(Γ (a ))2n! and rn = 0.

16.68.20.5 Jacobi Polynomials

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 deﬁnition is (−1, 1), the weight function w(x) = (1 −x)^a(1 + x)^b and, for the orthonormal case, the leading coeﬃcients are given by

(n + a + b + 1) ˜hn = -------n-------n 2 n!

where (n + a + b + 1)_n is the Pochhammer symbol. The parameters a and b should satisfy a > −1, b > −1. Also

˜ a+b+1 Γ (a-+-1)Γ (b-+-1) h0 = 2 Γ (a + b + 2) ˜hn = 2a+b+1 ----Γ-(n-+-a +-1)Γ (n-+-b-+-1)---- for n > 0 (2n + a + b + 1)Γ (n + a + b + 1)n! n(a − b) rn = ----------. 2n + a + b

The Legendre, Chebyshev and Gegenbauer polynomials are all, in fact, special cases of the Jacobi polynomials.

16.68.20.6 Laguerre Polynomials

The function call LaguerreP(n,x) will return the nth 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 deﬁnition is (0,∞), the weight function w(x) = e^−x and, for the orthonormal case, the leading coeﬃcients are given by k_n = (−1)ⁿ∕n!, h_n = 1 and r_n = −n².

16.68.20.7 Generalised Laguerre Polynomials

If used with three arguments LaguerreP(n,a,x) returns the nth 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 deﬁnition is (0,∞), the weight function w(x) = e^−xx^a and, for the orthonormal case, the leading coeﬃcients are given by k_n = (−1)ⁿ∕n!, h_n = Γ(n + a + 1)∕n! and r_n = −n(n + a). The parameter a should satisfy a > −1.

16.68.20.8 Hermite Polynomials

The function call HermiteP(n,x) will return the nth 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 deﬁnition is (−∞, +∞), the weight function w(x) = e^−x² and, for the orthonormal case, the leading coeﬃcients are given by k_n = 2ⁿ, h_n = √ --
π 2ⁿn! and r_n = 0.

16.68.21 Other Polynomials and Related Numbers

16.68.21.1 Fibonacci Polynomials

FibonacciP(n,x) returns the nth Fibonacci polynomial in the variable x. If n is an integer, it will be evaluated using the recursive deﬁnition:

F0(x) = 0; F1(x) = 1; Fn(x ) = xFn−1(x ) + Fn −2(x).

The recursion is, of course, optimised as a simple loop to avoid repeated computation of lower-order polynomials.

16.68.21.2 Euler Numbers and Polynomials

Euler numbers are computed by the unary operator Euler; the call Euler(n) returns the nth Euler number; all the odd Euler numbers are zero. The computation is derived directly from Pascal’s triangle of binomial coeﬃcients.

The Euler numbers and polynomials have the following generating functions:

t ∑∞ n xt ∑∞ n -2e----= Ent--, --e--- = En-(x)t- 1 + e2t n=0 n! 1 + et n=0 n!

respectively. Thus E₀ = 1 and E₁ = 0. Furthermore the numbers and polynomials are related by the equations:

( 1) ∑n (n ) Ek ( 1)n −k En = 2nEn -- , En (x) = -k- x − -- . 2 k=0 k 2 2

The Euler polynomials are evaluated for non-negative integer n by using the summation immediately above.

16.68.21.3 Bernoulli Numbers & Polynomials

The call Bernoulli(n) evaluates to the nth 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 oﬀ, 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:

t ∑∞ B tn text ∞∑ B (x)tn ------= -n---, ------= -n------ et − 1 n=0 n! et − 1 n=0 n!

respectively. Thus B₀ = 1 and B₁ = − 1
2

. Furthermore the numbers and polynomials are related by the equations:

( ) ∑n n n− k Bn = Bn (0 ), Bn (x) = Bkx . k=0 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 nth 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.

16.68.22 Acknowledgements

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.

16.68.23 Tables of Operators and Constants

Special Functions

Function	Operator

Si(z)	Si(z)
Si(z) − π∕2	s_i(z)
Ci(z)	Ci(z)
Shi(z)	Shi(z)
Chi(z)	Chi(z)
erf(z)	Erf(z)
1 − erf(z)	erfc(z)
Ei(z)	Ei(z)
Ei(log(z))	Li(z)
C(x)	Fresnel_C(x)
S(x)	Fresnel_S(x)

B(a,b)	Beta(a,b)
Γ(a)	Gamma(a)
normalized incomplete Beta I_x(a,b) =	iBeta(a,b,x)
normalized incomplete Gamma P(a,z) =	iGamma(a,z)
incomplete Gamma γ(a,z)	m_gamma(a,z)
(a)_k	Pochhammer(a,k)
ψ(z)	Psi(z)
ψ⁽ⁿ⁾(z)	Polygamma(n,z)

J_ν(z)	BesselJ(nu,z)
Y _ν(z)	BesselY(nu,z)
I_ν(z)	BesselI(nu,z)
K_ν(z)	BesselK(nu,z)
H_n⁽¹⁾(z)	Hankel1(n,z)
H_n⁽²⁾(z)	Hankel2(n,z)

More Special Functions

Function	Operator

Ai(z)	Airy_Ai(z)
Bi(z)	Airy_Bi(z)
Ai′(z)	Airy_Aiprime(z)
Bi′(z)	Airy_Biprime(z)
H_ν(z)	StruveH(nu,z)
L_ν(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 ₁F₁(a,b; z) or Φ(a,b; z)	KummerM(a,b,z)
U(a,b,z) or z^−a₂F₀(a,b; z) or Ψ(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 ζ(z)	zeta(z)
Lambert ω(z)	Lambert_W(z)
Li₂(z)	dilog(z)
Li_n(z)	Polylog(n,z)
Lerch’s transcendent Φ(z,s,a)	Lerch_Phi(z,s,a)

Function	Operator

Y _n^m(x,y,z,r2)	SolidHarmonicY(n,m,x,y,z,r2)
Y _n^m(𝜃,ϕ)	SphericalHarmonicY(n,m,theta,phi)

	ThreeJSymbol({j1,m1},{j2,m2},{j3,m3})

	Clebsch_Gordan({j1,m1},{j2,m2},{j3,m3})

	SixJSymbol({j1,j2,j3},{l1,l2,l3})

Elliptic Functions and Integrals

Function	Operator

sn(u,k)	jacobisn(u,k)
dn(u,k)	jacobidn(u,k)
cn(u,k)	jacobicn(u,k)
cd(u,k)	jacobicd(u,k)
sd(u,k)	jacobisd(u,k)
nd(u,k)	jacobind(u,k)
dc(u,k)	jacobidc(u,k)
nc(u,k)	jacobinc(u,k)
sc(u,k)	jacobisc(u,k)
ns(u,k)	jacobins(u,k)
ds(u,k)	jacobids(u,k)
cs(u,k)	jacobics(u,k)
am(u,k)	jacobiam(u,k)
Complete Integral (1st kind) K(k)	ellipticK(k)
K^′(k)	ellipticK!’(k)
Incomplete Integral (1st kind) F(ϕ,k)	ellipticF(phi,k)
Complete Integral (2nd kind) E(k)	ellipticE(k)
E^′(k)	ellipticE!’(k)
Legendre Incomplete Int (2nd kind) E(u,k)	ellipticE(u,k)
Jacobi Incomplete Int (2nd kind) ℰ(u,k)	JacobiE(u,k)
Jacobi’s Zeta Z(u,k)	jacobizeta(u,k)
𝜗₁(u,q)	elliptictheta1(u,q)
𝜗₂(u,q)	elliptictheta2(u,q)
𝜗₃(u,q)	elliptictheta3(u,q)
𝜗₄(u,q)	elliptictheta4(u,q)

Polynomial Functions

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^(α)(x)	GegenbauerP(n,alpha,x)
P_n^(α,β)(x)	JacobiP(n,alpha,beta,x)

Well-known Numbers and Reserved Constants

Function	Operator

	Binomial(n,m)

Fibonacci Numbers F_n	Fibonacci(n)
s_n^m	Stirling1(n,m)
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
Base of natural logarithms	e
Euler’s γ constant	Euler_gamma
Catalan’s constant	Catalan
Khinchin’s constant	Khinchin
Golden ratio	Golden_ratio

PrevTail

Front