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.

Groebner Bases for Modules

GSORT _ _ _ _ _ _ _ _ _ _ _ _ operator

syntax:

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

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

examples:

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

  y**2+2y*x+x**2

GSPLIT _ _ _ _ _ _ _ _ _ _ _ _ operator

syntax:

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.

examples:

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

  {y**2,2y*x+x**2}

GSPOLY _ _ _ _ _ _ _ _ _ _ _ _ operator

syntax:

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.

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.

syntax:

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

examples:

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.

syntax:

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

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

examples:

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.

syntax:

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.

examples:

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.

syntax:

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

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

examples:

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.

syntax:

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.

examples:

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.

syntax:

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

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

examples:

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.

syntax:

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

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

examples:

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.

syntax:

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.

syntax:

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.

syntax:

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.

syntax:

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 expressionsYou have tried to use A i n some way other than gamma5 in a high-energy physics expression.

Gamma5 not allowed unless vecdim is 4You 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 operatorA 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.

syntax:

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

examples:

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

syntax:

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_mintries 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 equations, where the variables are e quated to the coordinates of the result point.

examples:

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.

syntax:

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

examples:

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

syntax:

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.

examples:

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.

syntax:

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.

examples:

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.

syntax:

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

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

examples:

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.

syntax:

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.

examples:

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.

syntax:

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.

examples:

 
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.

syntax:

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.

examples:

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.

syntax:

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

examples:

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.

syntax:

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.

examples:

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.

syntax:

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.

syntax:

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.

examples:

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.

syntax:

roots(<p>)

where <p> is a univariate polynomial.

examples:

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)

examples:

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.

syntax:

examples:

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.

syntax:

examples:

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.

syntax:

examples:

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.

syntax:

examples:

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.