reduce-historical/r36/doc/ROOTS.TEX

\documentstyle[11pt,reduce]{article}
\title{The REDUCE Root Finding Package \\ Mod 1.94, 28 May 1993}
\date{}
\author {Stanley L. Kameny \\ E-mail: valley!stan@rand.org}
\begin{document}
\maketitle
\index{root finding} \index{ROOTS package}

\section{Introduction}

The root finding package is designed so that it can be used as an
independent package, or it can be integrated with and called by {\tt
SOLVE}. \index{SOLVE package ! with ROOTS package} This document describes
the package in its independent use.  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.

\section{Root Finding Strategies}

For all polynomials handled by the root finding package, strategies of
factoring are employed where possible to reduce the amount of required
work.  These include square-free factoring and separation of complex
polynomials into a product of a polynomial with real coefficients and one
with complex coefficients.  Whenever these succeed, the resulting smaller
polynomials are solved separately, except that the root accuracy takes
into account the possibility of close roots on different branches.  One
other strategy used where applicable is the powergcd method of reducing
the powers of the initial polynomial by a common factor, and deriving the
roots in two stages, as roots of the reduced power polynomial.  Again
here, the possibility of close roots on different branches is taken into
account.

\section{Top Level Functions}

The top level functions can be called either as symbolic operators from
algebraic mode, or they can be called directly from symbolic mode with
symbolic mode arguments.  Outputs are expressed in forms that print out
correctly in algebraic mode.


\subsection{Functions that refer to real roots only}

Three top level functions refer only to real roots.  Each of these
functions can receive 1, 2 or 3 arguments.

The first argument is the polynomial p, that can be complex and can
have multiple or zero roots.  If arg2 and arg3 are not present, all real
roots are found.  If the additional arguments are present, they restrict
the region of consideration.

\begin{itemize}
\item If arguments are (p,arg2) then
Arg2 must be POSITIVE or NEGATIVE.  If arg2=NEGATIVE then only
negative roots of p are included; if arg2=POSITIVE then only positive
roots of p are included. Zero roots are excluded.

\item If arguments are (p,arg2,arg3) then
\ttindex{EXCLUDE} \ttindex{POSITIVE} \ttindex{NEGATIVE} \ttindex{INFINITY}
Arg2 and Arg3 must be r (a real number) or  EXCLUDE r,  or a member of
the list POSITIVE, NEGATIVE, INFINITY, -INFINITY.  EXCLUDE r causes the
value r to be excluded from the region.  The order of the sequence
arg2, arg3 is unimportant.  Assuming that arg2 $\leq$ arg3 when both are
numeric, then

\begin{tabular}{l c l}
\{-INFINITY,INFINITY\} & is equivalent to & \{\} represents all roots; \\
\{arg2,NEGATIVE\} & represents & $-\infty < r < arg2$; \\
\{arg2,POSITIVE\} & represents & $arg2 < r < \infty$;
\end{tabular}

In each of the following, replacing an {\em arg} with EXCLUDE {\em arg}
converts the corresponding inclusive $\leq$ to the exclusive $<$

\begin{tabular}{l c l}
\{arg2,-INFINITY\} & represents & $-\infty < r \leq arg2$; \\
\{arg2,INFINITY\} & represents & $arg2 \leq r < \infty$; \\
\{arg2,arg3\} & represents & $arg2 \leq r \leq arg3$;
\end{tabular}

\item If zero is in the interval the zero root is included.
\end{itemize}

\begin{description}
\ttindex{REALROOTS} \index{Sturm Sequences}
\item[REALROOTS] This function finds the real roots of the polynomial p,
using the REALROOT package to isolate real roots by the method of Sturm
sequences, then polishing the root to the desired accuracy.  Precision
of computation is guaranteed to be sufficient to separate all real roots
in the specified region.  (cf. MULTIROOT for treatment of multiple
roots.)

\ttindex{ISOLATER}
\item[ISOLATER] This function produces a list of rational intervals, each
containing a single real root of the polynomial p, within the specified
region, but does not find the roots.

\ttindex{RLROOTNO}
\item[RLROOTNO] This function computes the number of real roots of p in
the specified region, but does not find the roots.
\end{description}

\subsection{Functions that return both real and complex roots}

\begin{description}
\ttindex{ROOTS}
\item[ROOTS p;] This 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.  The value returned by ROOTS is a
list of equations for all roots.  In addition, ROOTS stores separate lists
of real roots and complex roots in the global variables ROOTSREAL and
ROOTSCOMPLEX. \ttindex{ROOTSREAL} \ttindex{ROOTSCOMPLEX}

The order of root discovery by ROOTS is highly variable from system to
system, depending upon very subtle arithmetic differences during the
computation.  In order to make it easier to compare results obtained on
different computers, the output of ROOTS is sorted into a standard order:
a root with smaller real part precedes a root with larger real part; roots
with identical real parts are sorted so that larger imaginary part
precedes smaller imaginary part. (This is done so that for complex pairs,
the positive imaginary part is seen first.)

However, when a polynomial has been factored (by square-free factoring or
by separation into real and complex factors) then the root sorting is
applied to each factor separately.  This makes the final resulting order
less obvious.  However it is consistent from system to system.

\ttindex{ROOTS\_AT\_PREC}
\item[ROOTS\_AT\_PREC p;] Same as ROOTS except that roots values are
returned to a minimum of the number of decimal places equal to the current
system precision.

\ttindex{ROOT\_VAL}
\item[ROOT\_VAL p;] Same as ROOTS\_AT\_PREC, except that instead of
returning a list of equations for the roots, a list of the root value is
returned.  This is the function that SOLVE calls.

\ttindex{NEARESTROOT}
\item[NEARESTROOT(p,s);] This top level function uses an iterative method
to find the root to which the method converges given the initial starting
origin s, which can be complex.  If there are several roots in the
vicinity of s and s is not significantly closer to one root than it is to
all others, the convergence could arrive at a root that is not truly the
nearest root.  This function should therefore be used only when the user
is certain that there is only one root in the immediate vicinity of the
starting point s.

\ttindex{FIRSTROOT}
\item[FIRSTROOT p;] ROOTS is called, but only the first root determined by
ROOTS is computed.  Note that this is not in general the first root that
would be listed in ROOTS output, since the ROOTS outputs are sorted into
a canonical order.  Also, in some difficult root finding cases, the first
root computed might be incorrect.
\end{description}


\subsection{Other top level functions}

\begin{description}
\ttindex{GETROOT} \ttindex{ROOTS} \ttindex{REALROOTS} \ttindex{NEARESTROOTS}
\item[GETROOT(n,rr);] If rr has the form of the output of ROOTS, REALROOTS,
or NEARESTROOTS; GETROOT returns the rational, real, or complex value of
the root equation.  An error occurs if $n<1$ or $n>$ the number of roots in
rr.

\ttindex{MKPOLY}
\item[MKPOLY rr;] This function can be used to reconstruct a polynomial
whose root equation list is rr and whose denominator is 1.  Thus one can
verify that if $rr := ROOTS~p$, and $rr1 := ROOTS~MKPOLY~rr$, then
$rr1 = rr$. (This will be true if {\tt MULTIROOT} and {\tt RATROOT} are ON,
and {\tt ROUNDED} is off.)
However, $MKPOLY~rr - NUM~p = 0$ will be true if and only if all roots of p
have been computed exactly.

\end{description}

\subsection{Functions available for diagnostic or instructional use only}

\begin{description}
\ttindex{GFNEWT}
\item[GFNEWT(p,r,cpx);] This function will do a single pass through the
function GFNEWTON for polynomial p and root r.  If cpx=T, then any
complex part of the root will be kept, no matter how small.

\ttindex{GFROOT}
\item[GFROOT(p,r,cpx);] This function will do a single pass through the
function GFROOTFIND for polynomial p and root r.  If cpx=T, then any
complex part of the root will be kept, no matter how small.
\end{description}
\section{Switches Used in Input}

The input of polynomials in algebraic mode is sensitive to the switches
{\tt COMPLEX}, {\tt ROUNDED}, and {\tt ADJPREC}.  The correct choice of
input method is important since incorrect choices will result in
undesirable truncation or rounding of the input coefficients.

Truncation or rounding may occur if {\tt ROUNDED} is on and
one of the following is true:

\begin{enumerate}
\item a coefficient is entered in floating point form or rational form.
\item {\tt COMPLEX} is on and a coefficient is imaginary or complex.
\end{enumerate}

Therefore, to avoid undesirable truncation or rounding, then:

\begin{enumerate}
\item {\tt ROUNDED} should be off and input should be
in integer or rational form; or
\item {\tt ROUNDED} can be on if it is acceptable to truncate or round
input to the current value of system precision; or both {\tt ROUNDED} and
{\tt ADJPREC} can be on, in which case system precision will be adjusted
to accommodate the largest coefficient which is input; or \item if the
input contains complex coefficients with very different magnitude for the
real and imaginary parts, then all three switches {\tt ROUNDED}, {\tt
ADJPREC} and {\tt COMPLEX} must be on.
\end{enumerate}

\begin{description}
\item[integer and complex modes] (off {\tt ROUNDED}) any real
polynomial can be input using integer coefficients of any size; integer or
rational coefficients can be used to input any real or complex polynomial,
independent of the setting of the switch {\tt COMPLEX}.  These are the most
versatile input modes, since any real or complex polynomial can be input
exactly.

\item[modes rounded and complex-rounded] (on {\tt ROUNDED}) polynomials can be
input using
integer coefficients of any size.  Floating point coefficients will be
truncated or rounded, to a size dependent upon the system.  If complex
is on, real coefficients can be input to any precision using integer
form, but coefficients of imaginary parts of complex coefficients will
be rounded or truncated.
\end{description}

\section{Internal and Output Use of Switches}

The REDUCE arithmetic mode switches {\tt ROUNDED} and {\tt COMPLEX}
control the behavior of the root finding package.
These switches are returned in the same state in which they were set
initially, (barring catastrophic error).

\begin{description}
\ttindex{COMPLEX}
\item[COMPLEX] The root finding package controls the switch {\tt COMPLEX}
internally, turning the switch on if it is processing a complex
polynomial.
For a polynomial with real coefficients, the
\ttindex{NEARESTROOT}
starting point argument for NEARESTROOT can be given in algebraic mode in
complex form as rl + im * I and will be handled correctly, independent of
the setting of the switch {\tt COMPLEX.} Complex roots will be computed
and printed correctly regardless of the setting of the switch {\tt
COMPLEX}.  However, if {\tt COMPLEX} is off, the imaginary part will print
out ahead of the real part, while the reverse order will be obtained if
COMPLEX is on.

\ttindex{ROUNDED}
\item[ROUNDED] The
root finding package performs computations using the arithmetic mode that
is required at the time, which may be integer, Gaussian integer, rounded,
or complex rounded.  The switch {\tt BFTAG} is used internally to govern
the mode of computation and precision is adjusted whenever necessary.  The
initial position of switches {\tt ROUNDED} and {\tt COMPLEX} are ignored.
At output, these switches will emerge in their initial positions.
\end{description}

\section{Root Package Switches}

Note: switches {\tt AUTOMODE}, {\tt ISOROOT} and {\tt ACCROOT}, present in
earlier versions, have been eliminated.

\begin{description}
\ttindex{RATROOT}
\item[RATROOT] (Default OFF) If {\tt RATROOT} is on all root equations are
output in rational form.  Assuming that the mode is {\tt COMPLEX} (i.e.
{\tt ROUNDED} is off,) the root equations are
guaranteed to be able to be input into REDUCE without truncation or
rounding errors. (Cf. the function MKPOLY described above.)

\ttindex{MULTIROOT}
\item[MULTIROOT] (Default ON) Whenever the polynomial has complex
coefficients or has real coefficients and has multiple roots, as
\ttindex{SQFRF} determined by the Sturm function, the function {\tt SQFRF}
is called automatically to factor the polynomial into square-free factors.
If {\tt MULTIROOT} is on, the multiplicity of the roots will be indicated
in the output of ROOTS or REALROOTS by printing the root output
repeatedly, according to its multiplicity.  If {\tt MULTIROOT} is off,
each root will be printed once, and all roots should be normally be
distinct. (Two identical roots should not appear.  If the initial
precision of the computation or the accuracy of the output was
insufficient to separate two closely-spaced roots, the program attempts to
increase accuracy and/or precision if it detects equal roots.  If,
however, the initial accuracy specified was too low, and it was not
possible to separate the roots, the program will abort.)

\index{tracing ! ROOTS package}
\ttindex{TRROOT}
\item[TRROOT] (Default OFF) If switch {\tt TRROOT} is on, trace messages
are printed out during the course of root determination, to show the
progress of solution.

\ttindex{ROOTMSG} \item[ROOTMSG] (Default OFF) If switch
{\tt ROOTMSG} is on in addition to switch {\tt TRROOT,} additional
messages are printed out to aid in following the progress of Laguerre and
Newton complex iteration.  These messages are intended for debugging use
primarily.


\end{description}


\section{Operational Parameters and Parameter Setting.}

\begin{description}
\ttindex{ROOTACC\#}
\item[ROOTACC\#] (Default 6) This parameter can be set using the function
ROOTACC n; which causes {\tt ROOTACC\#} to be set to MAX(n,6).  If {\tt
ACCROOT} is on, roots will be determined to a minimum of {\tt ROOT\-ACC\#}
significant places. (If roots are closely spaced, a higher number of
significant places is computed where needed.)

\ttindex{system precision}
\item[system precision] The roots package, during its operation, will
change the value of system precision but will restore the original value
of system precision at termination except that the value of system
precision is increased if necessary to allow the full roots output to be
printed.

\ttindex{PRECISION}
\item[PRECISION n;] If the user sets system precision, using the command
PRECISION n; then the effect is to increase the system precision to n, and
to have the same effect on ROOTS as ROOTACC n; ie. roots will now be
printed with minimum accuracy n.  The original conditions can then be
restored by using the command PRECISION RESET; or PRECISION NIL;.

\ttindex{ROOTPREC}
\item[ROOTPREC n;] The roots package normally sets the computation mode and
precision automatically.  However, if ROOTPREC n; is
called and $n$ is greater than the initial system precision then all root
computation will be done initially using a minimum system precision n.
Automatic operation can be restored by input of ROOTPREC 0;.
\end{description}


\section{Avoiding truncation of polynomials on input}

The roots package will not internally truncate polynomials.  However, it
is possible that a polynomial can be truncated by input reading functions
of the embedding lisp system, particularly when input is given in floating
point (rounded) format.

To avoid any difficulties, input can be done in integer or Gaussian
integer format, or mixed, with integers or rationals used to represent
quantities of high precision. There are many examples of this in the
test package.  It is usually best to let the roots package
determine the precision needed to compute roots.

The number of digits that can be safely represented in floating point in
the lisp system are contained in the global variable {\tt !!NFPD}.
Similarly, the maximum number of significant figures in floating point
output are contained in the global variable {\tt !!FLIM}.  The roots
package computes these values, which are needed to control the logic of
the program. \ttindex{"!"!FLIM} \ttindex{"!"!NFPD}

The values of intermediate root iterations (that are printed when {\tt
TRROOT} is on) are given in bigfloat format even when the actual values
are computed in floating point.  This avoids intrusive rounding of root
printout.

\end{document}