**Tracing Groebner bases**

**MODULE**

Given a polynomial ring, e.g. R=Z[x,y,...] and an integer n>1. The vectors with n elements of R form a free MODULE under elementwise addition and multiplication with elements of R.

For a submodule given by a finite basis a Groebner basis
can be computed, and the facilities of the GROEBNER package
are available except the operators
groebnerf
and *groesolve*. The vectors are encoded using auxiliary
variables which represent the unit vectors in the module.
These are declared in the share variable
gmodule.

**GMODULE** _ _ _ _ _ _ _ _ _ _ _ _ **variable**

The vectors of a free
module over a polynomial ring R
are encoded as linear combinations with unit vectors of
M which are represented by auxiliary variables. These
must be collected in the variable *gmodule* before
any call to an operator of the Groebner package.

torder({x,y,v1,v2,v3})$ gmodule := {v1,v2,v3}$ g:=groebner({x^2*v1 + y*v2,x*y*v1 - v3,2y*v1 + y*v3});

compute the Groebner basis of the submodule

([x^2,y,0],[xy,0,-1],[0,2y,y])

The members of the list *gmodule* are automatically
appended to the end of the variable list, if they are not
yet members there. They take part in the actual term ordering.

**INDEX**

**Groebner Bases for Modules**

**GSORT** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

where <p> is a polynomial or a list of polynomials.

The polynomials are reordered and sorted corresponding to the current term order.

torder lex; gsort(x**2+2x*y+y**2,{y,x}); y**2+2y*x+x**2

**GSPLIT** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

where <p> is a polynomial or a list of polynomials.

The polynomial is reordered corresponding to the the current term order and then separated into leading term and reductum. Result is a list with the leading term as first and the reductum as second element.

torder lex; gsplit(x**2+2x*y+y**2,{y,x}); {y**2,2y*x+x**2}

**GSPOLY** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

where <p1> and <p2> are polynomials.

The *subtraction* polynomial of p1 and p2 is computed
corresponding to the method of the Buchberger algorithm for
computing *groebner bases*: p1 and p2 are multiplied
with terms such that when subtracting them the leading terms
cancel each other.

**INDEX**

**Computing with distributive polynomials**

**Groebner package**

**HEPHYS** _ _ _ _ _ _ _ _ _ _ _ _ **introduction**

The High-energy Physics package is historic for REDUCE, since REDUCE originated as a program to aid in computations with Dirac expressions. The commutation algebra of the gamma matrices is independent of their representation, and is a natural subject for symbolic mathematics. Dirac theory is applied to beta decay and the computation of cross-sections and scattering. The high-energy physics operators are available in the REDUCE main program, rather than as a module which must be loaded.

**.** _ _ _ **HE-DOT** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The . operator is used to denote the scalar product of two Lorentz four-vectors.

<vector> *.* <vector>

<vector> must be an identifier declared to be of type *vector* to h
ave
the scalar product definition. When applied to arguments that are not
vectors, the
cons operator is used,
whose symbol is also ``dot.''

vector aa,bb,cc; let aa.bb = 0; aa.bb; 0 aa.cc; AA.CC q := aa.cc; Q := AA.CC q; AA.CC

Since vectors are special high-energy physics entities that do not contain values, the . product will not return a true scalar product. You can assign a scalar identifier to the result of a . operation, or assign a . operation to have the value of the scalar you supply, as shown above. Note that the result of a . operation is a scalar, not a vector.

The metric tensor g(u,v) can be represented by *u.v*. If contraction
over the indices is required, *u* and *v* should be declared to
be of type
index.

The dot operator has the highest precedence of the infix operators, so expressions involving . and other operators have the scalar product evaluated first before other operations are done.

**EPS** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The *eps* operator denotes the completely antisymmetric tensor of
order 4 and its contraction with Lorentz four-vectors, as used in
high-energy physics calculations.

*eps*(<vector-expr>,<vector-expr>,<vector-expr>,
<vector-expr>)

<vector-expr> must be a valid vector expression, and may be an index.

vector g0,g1,g2,g3; eps(g1,g0,g2,g3); - EPS(G0,G1,G2,G3); eps(g1,g2,g0,g3); EPS(G0,G1,G2,G3); eps(g1,g2,g3,g1); 0

Vector identifiers are ordered alphabetically by REDUCE. When an o
dd number
of transpositions is required to restore the canonical order to the four
arguments of *eps*, the term is ordered and carries a minus sign. When an
even number of transpositions is required, the term is returned ordered and
positive. When one of the arguments is repeated, the value 0 is returned.
A contraction of the form
eps(_i j mu nu p_mu q_nu)
is represented by *eps(i,j,p,q)* when *i* and *j* have been
declared to be of type
index.

**G** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

*g* is an n-ary operator used to denote a product of gamma matrices
contracted with Lorentz four-vectors, in high-energy physics.

*g*(<identifier>,<vector-expr>
{,<vector-expr>}*)

<identifier> is a scalar identifier representing a fermion line identifier, <vector-expr> can be any valid vector expression, representing a vector or a gamma matrix.

vector aa,bb,cc; vector a; g(line1,aa,bb); AA.BB g(line2,aa,a); 0 g(id,aa,bb,cc); 0 g(li1,aa,bb) + k; AA.BB + K let aa.bb = m*k; g(ln1,aa)*g(ln1,bb); K*M g(ln1,aa)*g(ln2,bb); 0

The vector *A* is reserved in arguments of *g* to de
note the
special gamma matrix gamma_5. It must be declared to
be a vector before you use it.

Gamma matrix expressions are associated with fermion lines in a Feynman
diagram. If more than one line occurs in an expression, the gamma
matrices involved are separate (operating in independent spin space), as
shown in the last two example lines above. A product of gamma matrices
associated with a single line can be entered either as a single *g*
command with several vector arguments, or as products of separate *g*
commands each with a single argument.

While the product of vectors is not defined, the product, sum and difference of several gamma expressions are defined, as is the product of a gamma expression with a scalar. If an expression involving gamma matrices includes a scalar, the scalar is treated as if it were the product of itself with a unit 4 x 4 matrix.

Dirac expressions are evaluated by computing the trace of the expression using the commutation algebra of gamma matrices. The algorithms used are described in articles by J. S. R. Chisholm in <Il Nuovo Cimento X,> Vol. 30, p. 426, 1963, and J. Kahane, <Journal of Mathematical Physics>, Vol. 9, p. 1732, 1968. The trace is then divided by 4 to distinguish between the trace of a scalar and the trace of an expression that is the product of a scalar with a unit 4 x 4 matrix.

Trace calculations may be prevented over any line identifier by declaring it to be nospur. If it is later desired to evaluate th ese traces, the declaration can be undone with the spur declaration.

The notation of Bjorken and Drell, <Relativistic Quantum Mechanics,>
1964, is assumed in all operations involving gamma matrices. For an
example of the use of *g* in a calculation, see the <REDUCE
User's Manual>.

**INDEX** _ _ _ _ _ _ _ _ _ _ _ _ **declaration**

The declaration *index* flags a four-vector as an index for subsequent
high-energy physics calculations.

*index*<vector-id>{,<vector-id>}*

<vector-id> must have been declared of type *vector*.

vector aa,bb,cc; index uu; let aa.bb = 0; (aa.uu)*(bb.uu); 0 (aa.uu)*(cc.uu); AA.CC

Index variables are used to represent contraction over components of vectors when scalar products are taken by the . operator, as well as indicating contraction for the eps operator or metric tensor.

The special status of a vector as an index can be revoked with the declaration remind. The object remains a vector, however.

**MASS** _ _ _ _ _ _ _ _ _ _ _ _ **command**

The *mass* command associates a scalar variable as a mass with
the corresponding vector variable, in high-energy physics calculations.

*mass*<vector-var>*=*<scalar-var>
{,<vector-var>*=*<scalar-var>}*

<vector-var> can be a declared vector variable; *mass* will declare
it to be of type *vector* if it is not. This may override an existing
matrix variable by that name. <scalar-var> must be a scalar variable.

vector bb,cc; mass cc=m; mshell cc; cc.cc; 2 M

Once a mass has been attached to a vector with a *mass* dec
laration,
the
mshell declaration puts the associated partic
le ``on the mass
shell.'' Subsequent scalar (.) products of the vector with itself will be
replaced by the square of the mass expression.

**MSHELL** _ _ _ _ _ _ _ _ _ _ _ _ **command**

The *mshell* command puts particles on the mass shell in high-energy
physics calculations.

*mshell*<vector-var>{,<vector-var>}*

<vector-var> must have had a mass attached to it by a mass declaration.

vector v1,v2; mass v1=m,v2=q; mshell v1; v1.v1; 2 M v2.v2; V2.V2 mshell v2; v1.v1*v2.v2; 2 2 M *Q

Even though a mass is attached to a vector variable representing a
particle, the replacement does not take place until the *mshell*
declaration is given for that vector variable.

**NOSPUR** _ _ _ _ _ _ _ _ _ _ _ _ **declaration**

The *nospur* declaration prevents the trace calculation over the given
line identifiers in high-energy physics calculations.

*nospur*<line-id>{,<line-id>}*

<line-id> is a scalar identifier that will be used as a line identifier.

vector a1,b1,c1; g(line1,a1,b1)*g(line2,b1,c1); A1.B1*B1.C1 nospur line2; g(line1,a1,b1)*g(line2,b1,c1); A1.B1*G(LINE2,B1,C1)

Nospur declarations can be removed by making the declaration spur.

**REMIND** _ _ _ _ _ _ _ _ _ _ _ _ **declaration**

The *remind* declaration removes the special status of its arguments
as indices, which was set in the
index declaration, in
high-energy physics calculations.

*remind*<identifier>{,<identifier>}*

<identifier> must have been declared to be of type index.

**SPUR** _ _ _ _ _ _ _ _ _ _ _ _ **declaration**

The *spur* declaration removes the special exemption from trace
calculations that was declared by
nospur, in high-energy physics
calculations.

*spur*<line-id>{,<line-id>}*

<line-id> must be a line-identifier that has previously been declared
*nospur*.

**VECDIM** _ _ _ _ _ _ _ _ _ _ _ _ **command**

The command *vecdim* changes the vector dimension from 4 to an arbitrary
integer or symbol. Used in high-energy physics calculations.

*vecdim*<dimension>

<dimension> must be either an integer or a valid scalar identifier that does not have a floating-point value.

The
eps operator and the gamma_5
symbol (*A*) are not properly defined in anything except four
dimensions and will print an error message if you use them that way. The
other high-energy physics operators should work without problem.

**VECTOR** _ _ _ _ _ _ _ _ _ _ _ _ **declaration**

The *vector* declaration declares that its arguments are of type *vect
or*.

*vector*<identifier>{,<identifier>}*

<identifier> must be a valid REDUCE identifier. It may have already been used for a matrix, array, operator or scalar variable. After an identifier has been declared to be a vector, it may not be used as a scalar variable.

Vectors are special entities for high-energy physics calculations. You
cannot put values into their coordinates; they do not have coordinates.
They are legal arguments for the high-energy physics operators
eps,
g and *.* (dot). Vector variables are
used to represent gamma matrices and gamma matrices contracted with Lorentz
4-vectors, since there are no Dirac variables per se in the system.
Vectors do follow the usual vector rules for arithmetic operations:
*+* and *-* operate upon two or more vectors, producing a
vector; *** and */* cannot be used between vectors; the
scalar product is represented by the . operator; and the product of a
scalar and vector expression is well defined, and is a vector.

You can represent components of vectors by including representations of unit
vectors in your system. For instance, letting *E0* represent the unit
vector (1,0,0,0), the command

*V1.E0 := 0;*would set up the substitution of zero for the first componen
t of the vector
*V1*.

Identifiers that are declared by the *index* and *mass* declaratio
ns are
automatically declared to be vectors.

The following errors can occur in calculations using the high energy physics package:

*A represents only gamma5 in vector expressions*You have tried to use A i
n some way other than gamma5 in a
high-energy physics expression.

*Gamma5 not allowed unless vecdim is 4*You have used gamma_5 in a high-en
ergy physics
computation involving a vector dimension other than 4.

<ID> *has no mass*

One of the arguments to mshell has had no mass assigned to it, in high-energy physics calculations.

*Missing arguments for G operator*A line symbol is missing in a gamma mat
rix expression in high-energy physics
calculations.

*Unmatched index*<list>

The parser has found unmatched indices during the evaluation of a gamma matrix expression in high-energy physics calculations.

**High Energy Physics**

**NUMERIC PACKAGE** _ _ _ _ _ _ _ _ _ _ _ _ **introduction**

The numeric package supplies algorithms based on approximation techniques of numerical mathematics. The algorithms use the rounded mode arithmetic of REDUCE, including the variable precision feature which is exploited in some algorithms in an adaptive manner in order to reach the desired accuracy.

**INTERVAL** _ _ _ _ _ _ _ _ _ _ _ _ **type**

Intervals are generally coded as lower bound and
upper bound connected by the operator *..*, usually
associated to a variable in an
equation.

where <var> is a kernel and <low>, <high> are numbers or expression which evaluate to numbers with <low><=<high >.

means that the variable x is taken in the range from 2 .5 up to 3.5.

**NUMERIC ACCURACY**

The keyword parameters *accuracy=a* and *iterations=i*,
where *a*and *i* must be positive integer numbers, control the
iterative algorithms: the iteration is continued until
the local error is below 10**-a; if that is impossible
within *i* steps, the iteration is terminated with an
error message. The values reached so far are then returned
as the result.

**TRNUMERIC** _ _ _ _ _ _ _ _ _ _ _ _ **switch**

Normally the algorithms produce only a minimum of printed
output during their operation. In cases of an unsuccessful
or unexpected long operation a *trace of the iteration* can be
printed by setting *trnumeric* *on*.

**NUM_MIN** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The Fletcher Reeves version of the *steepest descent*
algorithms is used to find the *minimum* of a
function of one or more variables. The
function must have continuous partial derivatives with respect to all
variables. The starting point of the search can be
specified; if not, random values are taken instead.
The steepest descent algorithms in general find only local
minima.

or

*num_min*(exp, {
<var>[=<val>] [,<var>[=<val>] ...] }
[,accuracy=<a>] [,iterations=<i>])

where <exp> is a function expression, <var> are the variables in <exp> and <val> are the (optional) start values. For <a> and <i> see numeric accuracy.

*Num_min*tries to find the next local minimum along the descending
path starting at the given point. The result is a
list
with the minimum function value as first element followed by a list
of
equation*s*, where the variables are e
quated to the coordinates
of the result point.

num_min(sin(x)+x/5, x) {4.9489585606,{X=29.643767785}} num_min(sin(x)+x/5, x=0) { - 1.3342267466,{X= - 1.7721582671}}

**NUM_SOLVE** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

An adaptively damped Newton iteration is used to find an approximative root of a function (function vector) or the solution of an equation (equation system). The expressions must have continuous derivatives for all variables. A starting point for the iteration can be given. If not given random values are taken instead. When the number of forms is not equal to the number of variables, the Newton method cannot be applied. Then the minimum of the sum of absolute squares is located instead.

With complex on, solutions with imaginary parts ca n be found, if either the expression(s) or the starting point contain a nonzero imaginary part.

or

*num_solve*({<exp>,...,<exp>}, <var>[=<val>],...,
<var>[=<val>]
[,accuracy=<a>][,iterations=<i>])

or

*num_solve*({<exp>,...,<exp>}, {<var>[=<val>],...
,<var>[=<val>]}
[,accuracy=<a>][,iterations=<i>])

where <exp> are function expressions, <var> are the variables, <val> are optional start values. For <a> and <i> see numeric accuracy.

*num_solve*tries to find a zero/solution of the expression(s).
Result is a list of equations, where the variables are
equated to the coordinates of the result point.

The *Jacobian matrix* is stored as side effect the shared
variable *jacobian*.

num_solve({sin x=cos y, x + y = 1},{x=1,y=2}); {X= - 1.8561957251,Y=2.856195584} jacobian; [COS(X) SIN(Y)] [ ] [ 1 1 ]

**NUM_INT** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

For the numerical evaluation of univariate integrals over a finite interval the following strategy is used: If int finds a formal antiderivative which is bounded in the integration interval, this is evaluated and the end points and the difference is returned. Otherwise a Chebyshev fit is computed, starting with order 20, eventually up to order 80. If that is recognized as sufficiently convergent it is used for computing the integral by directly integrating the coefficient sequence. If none of these methods is successful, an adaptive multilevel quadrature algorithm is used.

For multivariate integrals only the adaptive quadrature is used.
This algorithm tolerates isolated singularities.
The value *iterations* here limits the number of
local interval intersection levels.
<a> is a measure for the relative total discretization
error (comparison of order 1 and order 2 approximations).

where <exp> is the function to be integrated, <var> are the integration variables, <l> are the lower bounds, <u> are the upper bounds.

Result is the value of the integral.

num_int(sin x,x=(0 .. 3.1415926)); 2.0000010334

**NUM_ODESOLVE** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The *Runge-Kutta* method of order 3 finds an approximate graph for
the solution of real *ODE initial value problem*.

or

*num_odesolve*({<exp>,<exp>,...},
{ <depvar>=<start>,<depvar>=<start>,...}
<indep>=(<from> .. <to>)
[,accuracy=<a>][,iterations=<i>])

where <depvar> and <start> specify the dependent variable(s) and the starting point value (vector), <indep>, <from> and <to> specify the independent variable and the integration interval (starting point and end point), <exp> are equations or expressions which contain the first derivative of the independent variable with respect to the dependent variable.

The ODEs are converted to an explicit form, which then is used for a Runge Kutta iteration over the given range. The number of steps is controlled by the value of <i> (default: 20). If the steps are too coarse to reach the desired accuracy in the neighborhood of the starting point, the number is increased automatically.

Result is a list of pairs, each representing a point of the approximate solution of the ODE problem.

depend(y,x); num_odesolve(df(y,x)=y,y=1,x=(0 .. 1), iterations=5); ,{0.2,1.2214},{0.4,1.49181796},{0.6,1.8221064563}, {0.8,2.2255208258},{1.0,2.7182511366}}

In most cases you must declare the dependency relation between the variables explicitly using depend; otherwise the formal derivative might be converted to zero.

The operator solve is used to convert the form into an explicit ODE. If that process fails or if it has no unique result, the evaluation is stopped with an error message.

**BOUNDS** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

Upper and lower bounds of a real valued function over an
interval or a rectangular multivariate domain
are computed
by the operator *bounds*. The algorithmic basis is the computation
with inequalities: starting from the interval(s) of the
variables, the bounds are propagated in the expression
using the rules for inequality computation. Some knowledge
about the behavior of special functions like ABS, SIN, COS, EXP, LOG,
fractional exponentials etc. is integrated and can be evaluated
if the operator *bounds* is called with rounded mode on
(otherwise only algebraic evaluation rules are available).

If *bounds* finds a singularity within an interval, the evaluation
is stopped with an error message indicating the problem part
of the expression.

or

*bounds*(<exp>,{<var>=(<l> .. <u>)
[,<var>=(<l> .. <u>) ...]})

where <exp> is the function to be investigated,
<var> are the variables of <exp>,
<l> and <u> specify the area as set of
interval*s*.

*bounds*computes upper and lower bounds for the expression in the
given area. An
interval is returned.

bounds(sin x,x=(1 .. 2)); -1 .. 1 on rounded; bounds(sin x,x=(1 .. 2)); 0.84147098481 .. 1 bounds(x**2+x,x=(-0.5 .. 0.5)); - 0.25 .. 0.75

**CHEBYSHEV FIT**

The operator family *Chebyshev_...* implements approximation
and evaluation of functions by the Chebyshev method.
Let *T(n,a,b,x)* be the Chebyshev polynomial of order *n*
transformed to the interval *(a,b)*.
Then a function *f(x)* can be
approximated in *(a,b)* by a series

for i := 0:n sum c(i)*T(i,a,b,x)

The operator *chebyshev_fit* computes this approximation an
d
returns a list, which has as first element the sum expressed
as a polynomial and as second element the sequence
of Chebyshev coefficients.
*Chebyshev_df* and *Chebyshev_int* transform a Chebyshev
coefficient list into the coefficients of the corresponding
derivative or integral respectively. For evaluating a Chebyshev
approximation at a given point in the basic interval the
operator *Chebyshev_eval* can be used.
*Chebyshev_eval* is based on a recurrence relation which is
in general more stable than a direct evaluation of the
complete polynomial.

*chebyshev_eval*(<coeffs>,<var>=(<lo> .. <hi>),
<var>=<pt>)

*chebyshev_df*(<coeffs>,<var>=(<lo> .. <hi>))

*chebyshev_int*(<coeffs>,<var>=(<lo> .. <hi>))

where <fcn> is an algebraic expression (the target function), <var> is the variable of <fcn>, <lo> and <hi> are numerical real values which describe an interval <lo> <<hi>, the integer <n> is the approximation order (set to 20 if missing), <pt> is a number in the interval and <coeffs> is a series of Chebyshev coefficients.

on rounded; w:=chebyshev_fit(sin x/x,x=(1 .. 3),5); 3 2 w := {0.03824*x - 0.2398*x + 0.06514*x + 0.9778, {0.8991,-0.4066,-0.005198,0.009464,-0.00009511}} chebyshev_eval(second w, x=(1 .. 3), x=2.1); 0.4111

**NUM_FIT** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The operator *num_fit* finds for a set of
points the linear combination of a given set of
functions (function basis) which approximates the
points best under the objective of the *least squares*
criterion (minimum of the sum of the squares of the deviation).
The solution is found as zero of the
gradient vector of the sum of squared errors.

where <vals> is a list of numeric values,
<var> is a variable used for the approximation,
<pts> is a list of coordinate values which correspond to
<var>,
<basis> is a set of functions varying in *var* which is used
for the approximation.

The result is a list containing as first element the function which approximates the given values, and as second element a list of coefficients which were used to build this function from the basis.

pts:=for i:=1 step 1 until 5 collect i$ vals:=for i:=1 step 1 until 5 collect for j:=1:i product j$ num_fit(vals,{1,x,x**2},x=pts); 2 {14.571428571*X - 61.428571429*X + 54.6,{54.6, - 61.428571429,14.571428571}}

**Numeric Package**

**ROOTS PACKAGE** _ _ _ _ _ _ _ _ _ _ _ _ **introduction**

The root finding package is designed so that it can be used to find some or all of the roots of univariate polynomials with real or complex coefficients, to the accuracy specified by the user.

Not all operators of *roots package* are described here. For using
the operators

*isolater*(intervals isolating real roots)

*rlrootno*(number of real roots in an interval)

*rootsat-prec*(roots at system precision)

*rootval*(result in equation form)

*firstroot*(computing only one root)

*getroot*(selecting roots from a collection)

please consult the full documentation of the package.

**MKPOLY** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

Given a roots list as returned by
roots,
the operator *mkpoly* constructs a
polynomial which has these numbers as roots.

*mkpoly*<rl>

where <rl> is a list with equations, which all have the same kernel on their left-hand sides and numbers as right-hand sides.

mkpoly{x=1,x=-2,x=i,x=-i}; x**4 + x**3 - x**2 + x - 2

Note that this polynomial is unique only up to a numeric factor.

**NEARESTROOT** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The operator *nearestroot* finds one root of a polynomial
with an iteration using a given starting point.

where <p> is a univariate polynomial and <pt> is a number.

nearestroot(x^2+2,2); {x=1.41421*i}

The minimal accuracy of the result values is controlled by rootacc.

**REALROOTS** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The operator *realroots* finds that real roots of a polynomial
to an accuracy that is sufficient to separate them and which is
a minimum of 6 decimal places.

*realroots*(<p>,<from>,<to>)

where <p> is a univariate polynomial.
The optional parameters <from> and <to> classify
an interval: if given, exactly the real roots in this
interval will be returned. <from> and <to>
can also take the values *infinity* or *-infinity*.
If omitted all real roots will be returned.
Result is a
list
of equations which represent the roots of the polynomial at the
given accuracy.

realroots(x^5-2); {x=1.1487} realroots(x^3-104*x^2+403*x-300,2,infinity); {x=3.0,x=100.0} realroots(x^3-104*x^2+403*x-300,-infinity,2); {x=1}

The minimal accuracy of the result values is controlled by rootacc.

**ROOTACC** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The operator *rootacc* allows you to set the accuracy
up to which the roots package computes its results.

*rootacc*(<n>)

Here <n> is an integer value. The internal accuracy of
the *roots* package is adjusted to a value of
*max(6,n)*. The default value is *6*.

**ROOTS** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The operator *roots*
is the main top level function of the roots package.
It will find all roots, real and complex, of the polynomial p
to an accuracy that is sufficient to separate them and which is
a minimum of 6 decimal places.

where <p> is a univariate polynomial. Result is a
list
of equations which represent the roots of the polynomial at the
given accuracy. In addition, *roots* stores
separate lists of real roots and complex roots in the global
variables
rootsreal and
rootscomplex.

roots(x^5-2); {x=-0.929316 + 0.675188*i, x=-0.929316 - 0.675188*i, x=0.354967 + 1.09248*i, x=0.354967 - 1.09248*i, x=1.1487}

The minimal accuracy of the result values is controlled by rootacc.

**ROOT\_VAL** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The operator *root_val* computes the roots of a
univariate polynomial at system precision
(or greater if required for root separation) and presents
its result as a list of numbers.

*roots*(<p>)

where <p> is a univariate polynomial.

root_val(x^5-2); {-0.929316490603 + 0.6751879524*i, -0.929316490603 - 0.6751879524*i, 0.354967313105 + 1.09247705578*i, 0.354967313105 - 1.09247705578*i, 1.148698355}

**ROOTSCOMPLEX** _ _ _ _ _ _ _ _ _ _ _ _ **variable**

When the operator
roots is called the complex
roots are collected in the global variable *rootscomplex*
as
list.

**ROOTSREAL** _ _ _ _ _ _ _ _ _ _ _ _ **variable**

When the operator
roots is called the real
roots are collected in the global variable *rootreal*
as
list.

**Roots Package**

**SPECIAL FUNCTION PACKAGE** _ _ _ _ _ _ _ _ _ _ _ _ **introduction**

The REDUCE *Special Function Package* supplies extended
algebraic and numeric support for a wide class of objects.
This package was released together with REDUCE 3.5 (October 1993)
for the first time, a major update is released with REDUCE 3.6.

The functions included in this package are in most cases (unless otherwise stated) defined and named like in the book by Abramowitz and Stegun: Handbook of Mathematical Functions, Dover Publications.

The aim is to collect as much information on the special functions and simplification capabilities as possible, i.e. algebraic simplifications and numeric (rounded mode) code, limits of the functions together with the definitions of the functions, which are in most cases a power series, a (definite) integral and/or a differential equation.

What can be found: Some famous constants, a variety of Bessel functions, special polynomials, the Gamma function, the (Riemann) Zeta function, Elliptic Functions, Elliptic Integrals, 3J symbols (Clebsch-Gordan coefficients) and integral functions.

What is missing: Mathieu functions, LerchPhi, etc.. The information about the special functions which solve certain differential equation is very limited. In several cases numerical approximation is restricted to real arguments or is missing completely.

The implementation of this package uses REDUCE rule sets to a large extent, which guarantees a high 'readability' of the functions definitions in the source file directory. It makes extensions to the special functions code easy in most cases too. To look at these rules it may be convenient to use the showrules operator e.g.

showrulesBesseli;

.

Some evaluations are improved if the special function package is loaded, e.g. some (infinite) sums and products leading to expressions including special functions are known in this case.

Note: The special function package has to be loaded explicitly by calling

load_package specfn;

The functions MeijerG and hypergeometric require additionally

load_package specfn2;

**CONSTANTS**

There are a few constants known to the special function package, namely

_ _ _ *Euler's constant *(which can be computed as -
Psi(1)) and

_ _ _ *Khinchin's constant* (which is defined in Khinchin's book
``Continued Fractions'') and

_ _ _ *Golden_Ratio* (which can be computed as (1 + sqrt 5)/2) and

_ _ _ *Catalan's constant* (which is known as an infinite sum of recipro
cal
powers)

on rounded; Euler_Gamma; 0.577215664902 Khinchin; 2.68545200107 Catalan 0.915965594177 Golden_Ratio 1.61803398875

**BERNOULLI** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The *bernoulli* operator returns the nth Bernoulli number.

bernoulli 20; - 174611 / 330 bernoulli 17; 0

All Bernoulli numbers with odd indices except for 1 are zero.

**BERNOULLIP** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The *BernoulliP* operator returns the nth Bernoulli Polynomial
evaluated at x.

BernoulliP(3,z); 2 z*(2*z - 3*z + 1)/2 BernoulliP(10,3); 338585 / 66

The value of the nth Bernoulli Polynomial at 0 is the nth Bernoull i number.

**EULER** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The *EULER* operator returns the nth Euler number.

Euler 20; 370371188237525 Euler 0; 1

The *Euler* numbers are evaluated by a recursive algorithm
which
makes it hard to compute Euler numbers above say 200.

Euler numbers appear in the coefficients of the power series representation of 1/cos(z).

**EULERP** _ _ _ _ _ _ _ _ _ _ _ _ **operator**

The *EulerP* operator returns the nth Euler Polynomial.

EulerP(2,xx); xx*(xx - 1) EulerP(10,3); 2046

The Euler numbers are the values of the Euler Polynomials at 1/2 multiplied by 2**n.