REDUCE

16.67 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 first one is more general in nature, whereas the second is devoted to special special functions. Documentation for the first package can be found in the file specfn.tex in the “doc” directory, and examples in specfn.tst and specfmor.tst in the examples directory.

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

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

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 file specplot.tst in the subfolder plot of the packages folder. These examples will reproduce a number of well-known pictures from [AS72].

16.67.1 Simplification and Approximation

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 later package contains some information which may be useful and directly applicable for other packages, e.g.:

sum(1/k^s,k,1,infinity); % will be evaluated 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.

16.67.2 Constants

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.

16.67.3 Bernoulli Numbers and Euler Numbers

The unary operator Bernoulli provides notation and computation for Bernoulli numbers. Bernoulli(n) evaluates to the nth Bernoulli number; all of the odd Bernoulli numbers, except Bernoulli(1), are zero.

The algorithms 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.

Euler numbers are computed by the unary operator Euler, which return the nth Euler number. The computation is derived directly from Pascal’s triangle of binomial coefficients.

16.67.4 Fibonacci Numbers and Fibonacci Polynomials

The unary operator Fibonacci provides notation and computation for Fibonacci numbers. Fibonacci(n) evaluates to the nth Fibonacci number. If n is a positive or negative integer, it will be evaluated following the definition:

F0 = 0;F1 =  1;Fn = Fn −1 + Fn−2

Fibonacci Polynomials are computed by the binary operator FibonacciP. FibonacciP(n,x) returns the nth Fibonaccip polynomial in the variable x. If n is a positive or negative integer, it will be evaluated following the definition:

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

16.67.5 Stirling Numbers

Stirling numbers of the first and second kind are computed by the binary operators Stirling1 and Stirling2 using explicit formulae.

16.67.6 The Γ Function, and Related Functions

16.67.6.1 The Γ Function

This is represented by the unary operator Gamma.

Initial transformations applied with ROUNDED off are: Γ(n) for integral n is computed, Γ(n + 12) 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 functions.

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

16.67.6.2 The Pochhammer Notation

The Pochhammer notation (a)k is supported by the binary operator Pochhammer. With ROUNDED off, this expression is evaluated numerically if a and k are both integral, and otherwise may be simplified where appropriate. The simplification rules are based upon algorithms supplied by Wolfram Koepf [Koe92].

16.67.6.3 The Digamma Function, ψ

This is represented by the unary operator PSI.

Initial transformations for ψ are applied on a similar basis to those for Γ; where possible, ψ(x) is rewritten in terms of ψ(1) and ψ(1
2), 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 scaler.

Relations for the derivative and integral of ψ are included.

16.67.6.4 The Polygamma Functions, ψ(n)

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

Initial manipulations on ψ(n) are few; where the second argument is 1 or 32, the expression is rewritten to one involving the Riemann ζ function, and when the first is zero it is rewritten to ψ; poles are also handled.

Numerical evaluation is avaialbe 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 nth derivative of cot().

Simple relations for derivatives and integrals are provided.

16.67.6.5 The Riemann ζ Function

This is represented by the unary operator Zeta.

With ROUNDED off, ζ(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 ψn1(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 ζ.

16.67.7 Bessel Functions

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 following initial transformations are performed:

Also, if the COMPLEX switch is on and ROUNDED is off, expressions in Hankel functions are rewritten in terms of Bessels.

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 Bessels 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.

16.67.8 Hypergeometric and Other Functions

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

16.67.9 Integral Functions

The SPECFN package includes manipulation and a limited numerical evaluation for some Integral functions, namely

erf, erfc, Si, Shi, si, Ci, Chi, Ei, li, Fresnel_C, and Fresnel_S.

The definitions from integral, the derviatives and some limits are known together with some simple properties such as symmetry conditions.

The numerical approximation 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.

li(z) is simplified towards Ei(ln(z)) .

16.67.10 Airy Functions

Support is provided for the Airy Functions Ai and Bi and for the Airyprime Functions Aiprime and Biprime. The relevant operators are respectively Airy_Ai, Airy_Bi, Airy_Aiprime, and Airy_Biprime, which are all unary operators with one argument.

The following cases can be performed:

In order to activate the Airy Function to Bessel Rules one should type:
let Airy2Bessel_rules;. As a result the Airy_Ai function, for example will be calculated using the formula :-

         1√ -[                 ]             2 2
Ai(z ) = 3- z I− 1∕3(ζ) − I1∕3(ζ)  , where ζ = 3z3

Note:- In order to obtain satisfactory approximations to results 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 approach, rather than the series. This value is calculated by the program and depends on the given precision.

There are no rules for the integration of Airy Functions.

16.67.11 Polynomial Functions

Two groups are defined, some well-known orthogonal Polynomials (Hermite, Jacobi, Legendre, Laguerre, Chebyshev, Gegenbauer) and Euler and Bernoulli Polynomials. The names of the REDUCE operator are built by adding a P to the name of the polynomials, e.g. EulerP implements the Euler polynomials. Most definitions are equivalent to [AS72], except for the ternary (associated) Legendre Polynomials.

P(n,m,x) = (-1)^m *(1-x^2)^(m/2)*df(legendreP (n,x),x,m)

16.67.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 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 Laplace equation:-

       df(P,x,2) + df(P,y,2) + df(P,z,2) = 0.

There are 2*n+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 exp(i*k*u) with k neq 0,

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

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

Let r2 = x2 + y2 + z2 and r = sqrt(r2).

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 n >= 0; n <= m =< n form a complete orthogonal basis on it, i.e. < n,m|n,m> = Kronecker_delta(n,n’) * Kronecker_delta(m,m’) using < ...|... > 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 ortho-normal 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),...,(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);

16.67.13 Jacobi’s 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:

16.67.13.1 Jacobi Elliptic Functions

The following Jacobi functions are available:-

These differ 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).

Extended rule lists are provided for differentiation 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 iKand finally Jacobi’s transformation for a purely imaginary first 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.67.13.2 Jacobi Amplitude Function

The amplitude of u can be evaluated using the jacobiam(u,k) command.

16.67.13.3 Arithmetic Geometric Mean

A procedure to evaluate the AGM of initial values a0,b0,c0  exists as
AGM_function(a0,b0,c0     ) and will return
{N,AGM,{aN,,a0},{bN,,b0},{cN,,c0}}, where N is the number of steps to compute the AGM to the desired acuracy.

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

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.67.13.4 Descending Landen Transformation

The procedure to evaluate the Descending Landen Transformation of phi and alpha 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(ϕ0,α0) and will return
{{ϕ0,n},{α0,n}}.

16.67.14 Elliptic Integrals

The following functions have been implemented:

These again differ 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 first 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 avoid any possible confusion, following Lawden [Law89], the notation D(ϕ,k) will be used for the Legendre form.

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 simplifications for special values of the arguments.

16.67.14.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.67.14.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
K(k ) =     (1 − k2 sin2 𝜃)−1∕2d𝜃.
         0

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

16.67.14.3 Elliptic K

The Elliptic Kfunction can be used as EllipticK!(k) and will return the value K(k) = K(√ -------
  1 − k2). iK(k) is the other quarter-period of the Jacobi elliptic functions.

16.67.14.4 Elliptic D

The Elliptic D function takes two arguments. It is called as EllipticD(phi,k) and will return the value of Legendre’s form of the Incomplete Elliptic Integral of the Second Kind:

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

16.67.14.5 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 Jacobi’s form of the Incomplete Elliptic Integral of the Second Kind:

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

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

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

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

                     ∫ K(k)
E (k) = E(K (k),k) =       dn2(v, k)dv
                      0

or equivalently

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

16.67.14.6 Elliptic E

The Elliptic Efunction can be used as EllipticE!(k) and will return the value E(k) = E(√ -------
  1 − k2).

16.67.14.7 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.67.15 Elliptic Theta Functions

These theta functions differ from those originally defined by Lisa Temme in a number of respects. Firstly four separate functions of two arguments are defined

rather than a single function with three arguments (with the first 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 q14 in their defining 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 q14 is interpreted as exp(iπτ∕4) rather than the principal value of q14. Thus τ, 2 + τ, 4 + τ and 6 + τ produce four different 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 first argument into the base period parallelogram with vertices

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

Together with the relation 𝜗1(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 q14). 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 suffer 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 effect 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 qfor which |q′|≈|q|α2 where α = π∕| arg(q)|≥ 2. Thus the worst case occurs for values of the nome q near to ±i where α 2.

16.67.15.1 Some Numerical Utility Functions

Five utility functions are provided:

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

The second pair return the quarter-periods K and Krespectively 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.67.16 Lambert’s W function

Lambert’s W function is the inverse of the function w ew. Therefore it is an important contribution for the solve package. The function is studied extensively in [HCGJ92]. The current implementation will compute the principal branch in ROUNDED mode only.

16.67.17 3j symbols and Clebsch-Gordan Coefficients

The operators ThreeJSymbol, Clebsch_Gordan are defined like in [LB68] or [Edm57]. ThreeJSymbol expects as arguments three lists of values {ji,mi}, e.g.

ThreeJSymbol({J+1,M},{J,-M},{1,0});  
Clebsch_Gordan({2,0},{2,0},{2,0});

16.67.18 6j symbols

The operator SixJSymbol is defined like in [LB68] or [Edm57]. SixJSymbol expects two lists of values {j1,j2,j3} and {l1,l2,l3} as arguments, e.g.

SixJSymbol({7,6,3},{2,4,6});

In the current implementation of the SixJSymbol, there is only a 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.67.19 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.67.20 Table of Operators and Constants

FunctionOperator
  
Jν(z)BesselJ(nu,z)
Y ν(z)BesselY(nu,z)
Iν(z)BesselI(nu,z)
Kν(z)BesselK(nu,z)
Hν(1)(z)Hankel1(n,z)
Hν(2)(z)Hankel2(n,z)
Hν(z)StruveH(nu,z)
Lν(z)StruveL(n,z)
sa,b(z)Lommel1(a,b,z)
Sa,b(z)Lommel2(a,b,z)
Ai(z)Airy_Ai(z)
Bi(z)Airy_Bi(z)
Ai(z)Airy_Aiprime(z)
Bi(z)Airy_Biprime(z)
M(a,b,z) or 1F1(a,b; z) or Φ(a,b; z)KummerM(a,b,z)
U(a,b,z) or za 2F0(a,b; z) or Ψ(a,b; z)KummerU(a,b,z)
Mκ,μ(z)WhittakerM(kappa,mu,z)
Wκ,μ(z)WhittakerW(kappa,mu,z)
  
Fibonacci Numbers FnFibonacci(n)
Fibonacci Polynomials Fn(x)FibonacciP(n)
Bn(x)BernoulliP(n,x)
En(x)EulerP(n,x)
Cn(α)(x)GegenbauerP(n,alpha,x)
Hn(x)HermiteP(n,x)
Ln(x)LaguerreP(n,x)
Ln(m)(x)LaguerreP(n,m,x)
Pn(x)LegendreP(n,x)
Pn(m)(x)LegendreP(n,m,x)
FunctionOperator
  
Pn(α,β)(x)JacobiP(n,alpha,beta,x)
Un(x)ChebyshevU(n,x)
Tn(x)ChebyshevT(n,x)
Y nm(x,y,z,r2)SolidHarmonicY(n,m,x,y,z,r2)
Y nm(𝜃,ϕ)SphericalHarmonicY(n,m,theta,phi)
  
(              )
   j1   j2  j3
  m1   m2   m3ThreeJSymbol({j1,m1},{j2,m2},{j3,m3})
  
(j1m1j2m2 |j1j2j3 − m3 )Clebsch_Gordan({j1,m1},{j2,m2},{j3,m3})
  
{  j    j   j  }
    1    2   3
  m1   m2   m3SixJSymbol({j1,j2,j3},{l1,l2,l3})
  
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)
F(ϕ,k)ellipticF(phi,k)
K(k)ellipticK(k)
K(k)ellipticK!’(k)
D(ϕ,k)ellipticD(phi,k)
E(u,k)orE(k)ellipticE(u,k) or ellipticE(k)
E(k)ellipticE!’(k)
Z(u,k)jacobizeta(u,k)
𝜗1(u,q)elliptictheta1(u,q)
𝜗2(u,q)elliptictheta2(u,q)
𝜗3(u,q)elliptictheta3(u,q)
𝜗4(u,q)elliptictheta4(u,q)
Lambert ω(z)Lambert_W(z)
  
ConstantREDUCE name
Euler’s γ constantEuler_gamma
Catalan’s constantCatalan
Khinchin’s constant Khinchin
Golden ratioGolden_ratio

  
FunctionOperator
  
( n )

  mBinomial(n,m)
  
Motzkin(n)Motzkin(n)
Bernoulli(n) or BnBernoulli(n)
Euler(n) or EnEuler(n)
Sn(m)Stirling1(n,m)
Sn(m)Stirling2(n,m)
B(z,w)Beta(z,w)
Γ(z)Gamma(z)
normalized incomplete Beta Ix(a,b) = Bx(a, b)
B-(a,b)-iBeta(a,b,x)
normalized incomplete Gamma P(a,z) = γ(a,z)
-------
 Γ (a)iGamma(a,z)
incomplete Gamma γ(a,z)m_gamma(a,z)
(a)kPochhammer(a,k)
ψ(z)Psi(z)
ψ(n)(z)Polygamma(n,z)
Riemann’s ζ(z)Zeta(z)
  
Si(z)Si(z)
si(z)s_i(z)
Ci(z)Ci(z)
Shi(z)Shi(z)
Chi(z)Chi(z)
erf(z)erf(z)
erfc(z)erfc(z)
Ei(z)Ei(z)
li(z)li(z)
C(x)Fresnel_C(x)
S(x)Fresnel_S(x)
  
dilog(z)dilog(z)
Lin(z)Polylog(n,z)
Lerch’s transcendent Φ(z,s,a)Lerch_Phi(z,s,a)