Artifact 94266dad484ea1a26f4831e2672277c8a52ace3a3eb2d857624bfda80d9f747c:
- File
r36/doc/RANDPOLY.TEX
— part of check-in
[152fb3bdbb]
at
2011-10-17 17:58:33
on branch master
— svn:eol-style, svn:executable and line endings for files
in historical/r36 treegit-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1480 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: schoepf@users.sourceforge.net, size: 21045) [annotate] [blame] [check-ins using] [more...]
\documentstyle[11pt]{article} \title{RANDPOLY: A Random Polynomial Generator} \author{Francis J. Wright \\ School of Mathematical Sciences \\ Queen Mary and Westfield College \\ University of London \\ Mile End Road, London E1 4NS, UK. \\ Email: {\tt F.J.Wright@QMW.ac.uk}} \date{14 July 1994} \begin{document} \maketitle \begin{abstract} This package is based on a port of the Maple random polynomial generator together with some support facilities for the generation of random numbers and anonymous procedures. \end{abstract} \section{Introduction} The operator {\tt randpoly} is based on a port of the Maple random polynomial generator. In fact, although by default it generates a univariate or multivariate polynomial, in its most general form it generates a sum of products of arbitrary integer powers of the variables multiplied by arbitrary coefficient expressions, in which the variable powers and coefficient expressions are the results of calling user-supplied functions (with no arguments). Moreover, the ``variables'' can be arbitrary expressions, which are composed with the underlying polynomial-like function. The user interface, code structure and algorithms used are essentially identical to those in the Maple version. The package also provides an analogue of the Maple {\tt rand} random-number-generator generator, primarily for use by {\tt randpoly}. There are principally two reasons for translating these facilities rather than designing comparable facilites anew: (1) the Maple design seems satisfactory and has already been ``proven'' within Maple, so there is no good reason to repeat the design effort; (2) the main use for these facilities is in testing the performance of other algebraic code, and there is an advantage in having essentially the same test data generator implemented in both Maple and REDUCE\@. Moreover, it is interesting to see the extent to which a facility can be translated without change between two systems. (This aspect will be described elsewhere.) Sections \ref{sec:Basic} and \ref{sec:Advanced} describe respectively basic and more advanced use of {\tt randpoly}; \S\ref{sec:Subsidiary} describes subsidiary functions provided to support advanced use of {\tt randpoly}; \S\ref{sec:Examples} gives examples; an appendix gives some details of the only non-trivial algorithm, that used to compute random sparse polynomials. Additional examples of the use of {\tt randpoly} are given in the test and demonstration file {\tt randpoly.tst}. \section{Basic use of {\tt randpoly}} \label{sec:Basic} The operator {\tt randpoly} requires at least one argument corresponding to the polynomial variable or variables, which must be either a single expression or a list of expressions.% \footnote{If it is a single expression then the univariate code is invoked; if it is a list then the multivariate code is invoked, and in the special case of a list of one element the multivariate code is invoked to generate a univariate polynomial, but the result should be indistinguishable from that resulting from specifying a single expression not in a list.} % In effect, {\tt randpoly} replaces each input expression by an internal variable and then substitutes the input expression for the internal variable in the generated polynomial (and by default expands the result as usual), although in fact if the input expression is a REDUCE kernel then it is used directly. The rest of this document uses the term ``variable'' to refer to a general input expression or the internal variable used to represent it, and all references to the polynomial structure, such as its degree, are with respect to these internal variables. The actual degree of a generated polynomial might be different from its degree in the internal variables. By default, the polynomial generated has degree 5 and contains 6 terms. Therefore, if it is univariate it is dense whereas if it is multivariate it is sparse. \subsection{Optional arguments} Other arguments can optionally be specified, in any order, after the first compulsory variable argument. All arguments receive full algebraic evaluation, subject to the current switch settings etc. The arguments are processed in the order given, so that if more than one argument relates to the same property then the last one specified takes effect. Optional arguments are either keywords or equations with keywords on the left. In general, the polynomial is sparse by default, unless the keyword {\tt dense} is specified as an optional argument. (The keyword {\tt sparse} is also accepted, but is the default.) The default degree can be changed by specifying an optional argument of the form \begin{center} {\tt degree = {\it natural number}}. \end{center} In the multivariate case this is the total degree, i.e.\ the sum of the degrees with respect to the individual variables. The keywords {\tt deg} and {\tt maxdeg} can also be used in place of {\tt degree}. More complicated monomial degree bounds can be constructed by using the coefficient function described below to return a monomial or polynomial coefficient expression. Moreover, {\tt randpoly} respects internally the REDUCE ``asymptotic'' commands {\tt let}, {\tt weight} etc.\ described in \S10.4 of the \REDUCE 3.6 manual, which can be used to exercise additional control over the polynomial generated. In the sparse case (only), the default maximum number of terms generated can be changed by specifying an optional argument of the form \begin{center} {\tt terms = {\it natural number}}. \end{center} The actual number of terms generated will be the minimum of the value of {\tt terms} and the number of terms in a dense polynomial of the specified degree, number of variables, etc. \section{Advanced use of {\tt randpoly}} \label{sec:Advanced} The default order (or minimum or trailing degree) can be changed by specifying an optional argument of the form \begin{center} {\tt ord = {\it natural number}}. \end{center} The keyword is {\tt ord} rather than {\tt order} because {\tt order} is a reserved command name in REDUCE\@. The keyword {\tt mindeg} can also be used in place of {\tt ord}. In the multivariate case this is the total degree, i.e.\ the sum of the degrees with respect to the individual variables. The order normally defaults to 0. However, the input expressions to {\tt randpoly} can also be equations, in which case the order defaults to 1 rather than 0. Input equations are converted to the difference of their two sides before being substituted into the generated polynomial. The purpose of this facility is to easily generate polynomials with a specified zero -- for example \begin{center}\tt randpoly(x = a); \end{center} generates a polynomial that is guaranteed to vanish at $x = a$, but is otherwise random. Order specification and equation input are extensions of the current Maple version of {\tt randpoly}. The operator {\tt randpoly} accepts two further optional arguments in the form of equations with the keywords {\tt coeffs} and {\tt expons} on the left. The right sides of each of these equations must evaluate to objects that can be applied as functions of no variables. These functions should be normal algebraic procedures (or something equivalent); the {\tt coeffs} procedure may return any algebraic expression, but the {\tt expons} procedure must return an integer (otherwise {\tt randpoly} reports an error). The values returned by the functions should normally be random, because it is the randomness of the coefficients and, in the sparse case, of the exponents that makes the constructed polynomial random. A convenient special case is to use the function {\tt rand} on the right of one or both of these equations; when called with a single argument {\tt rand} returns an anonymous function of no variables that generates a random integer. The single argument of {\tt rand} should normally be an integer range in the form $a~..~b$, where $a$, $b$ are integers such that $a < b$. The spaces around (or at least before) the infix operator ``..'' are necessary in some cases in REDUCE and generally recommended. For example, the {\tt expons} argument might take the form \begin{center}\tt expons = rand(0~..~n) \end{center} where {\tt n} will be the maximum degree with respect to each variable {\em independently}. In the case of {\tt coeffs} the lower limit will often be the negative of the upper limit to give a balanced coefficient range, so that the {\tt coeffs} argument might take the form \begin{center}\tt coeffs = rand(-n~..~n) \end{center} which will generate random integer coefficients in the range $[-n,n]$. \section{Subsidiary functions: rand, proc, random} \label{sec:Subsidiary} \subsection{Rand: a random-number-generator generator} The first argument of {\tt rand} must be either an integer range in the form $a~..~b$, where $a$, $b$ are integers such that $a < b$, or a positive integer $n$ which is equivalent to the range $0~..~n-1$. The operator {\tt rand} constructs a function of no arguments that calls the REDUCE random number generator function {\tt random} to return a random integer in the range specified; in the case that the first argument of {\tt rand} is a single positive integer $n$ the function constructed just calls {\tt random($n$)}, otherwise the call of {\tt random} is scaled and shifted. As an additional convenience, if {\tt rand} is called with a second argument that is an identifier then the call of {\tt rand} acts exactly like a procedure definition with the identifier as the procedure name. The procedure generated can then be called with an empty argument list by the algebraic processor. [Note that {\tt rand()} with no argument is an error in REDUCE and does not return directly a random number in a default range as it does in Maple -- use instead the REDUCE function {\tt random} (see below).] \subsection{Proc: an anonymous procedure generator} The operator {\tt proc} provides a generalization of {\tt rand}, and is primarily intended to be used with expressions involving the {\tt random} function (see below). Essentially, it provides a mechanism to prevent functions such as {\tt random} being evaluated when the arguments to {\tt randpoly} are evaluated, which is too early. {\tt Proc} accepts a single argument which is converted into the body of an anonymous procedure, which is returned as the value of {\tt proc}. (If a named procedure is required then the normal REDUCE {\tt procedure} statement should be used instead.) Examples are given in the following sections, and in the file {\tt randpoly.tst}. \subsection{Random: a generalized interface} As an additional convenience, this package extends the interface to the standard REDUCE {\tt random} function so that it will directly accept either a natural number or an integer range as its argument, exactly as for the first argument of {\tt rand}. Hence effectively \begin{center}\tt rand(X) = proc random(X) \end{center} although {\tt rand} is marginally more efficient. However, {\tt proc} and the generalized {\tt random} interface allow expressions such as the following anonymous random fraction generator to be easily constructed: \begin{center}\tt proc(random(-99~..~99)/random(1~..~99)) \end{center} \subsection{Further support for procs} {\tt Rand} is a special case of {\tt proc}, and (for either) if the switch {\tt comp} is {\tt on} (and the compiler is available) then the generated procedure body is compiled. {\tt Rand} with a single argument and {\tt proc} both return as their values anonymous procedures, which if they are not compiled are Lisp lambda expressions. However, if compilation is in effect then they return only an identifier that has no external significance% \footnote{It is not interned on the oblist.} % but which can be applied as a function in the same way as a lambda expression. It is primarily intended that such ``proc expressions'' will be used immediately as input to {\tt randpoly}. The algebraic processor is not intended to handle lambda expressions. However, they can be output or assigned to variables in algebraic mode, although the output form looks a little strange and is probably best not displayed. But beware that lambda expressions cannot be evaluated by the algebraic processor (at least, not without declaring some internal Lisp functions to be algebraic operators). Therefore, for testing purposes or curious users, this package provides the operators {\tt showproc} and {\tt evalproc} respectively to display and evaluate ``proc expressions'' output by {\tt rand} or {\tt proc} (or in fact any lambda expression), in the case of {\tt showproc} provided they are not compiled. \section{Examples} \label{sec:Examples} The file {\tt randpoly.tst} gives a set of test and demonstration examples. The following additional examples were taken from the Maple {\tt randpoly} help file and converted to REDUCE syntax by replacing [~] by \{~\} and making the other changes shown explicitly: \begin{verbatim} randpoly(x); 5 4 3 2 - 54*x - 92*x - 30*x + 73*x - 69*x - 67 randpoly({x, y}, terms = 20); 5 4 4 3 2 3 3 31*x - 17*x *y - 48*x - 15*x *y + 80*x *y + 92*x 2 3 2 2 4 3 2 + 86*x *y + 2*x *y - 44*x + 83*x*y + 85*x*y + 55*x*y 5 4 3 2 - 27*x*y + 33*x - 98*y + 51*y - 2*y + 70*y - 60*y - 10 randpoly({x, sin(x), cos(x)}); 4 3 3 sin(x)*( - 4*cos(x) - 85*cos(x) *x + 50*sin(x) 2 - 20*sin(x) *x + 76*sin(x)*x + 96*sin(x)) % randpoly(z, expons = rand(-5..5)); % Maple % A generalized random "polynomial"! % Note that spaces are needed around .. in REDUCE. on div; off allfac; randpoly(z, expons = rand(-5 .. 5)); 4 3 -3 -4 -5 - 39*z + 14*z - 77*z - 37*z - 8*z off div; on allfac; % randpoly([x], coeffs = proc() randpoly(y) end); % Maple randpoly({x}, coeffs = proc randpoly(y)); 5 5 5 4 5 3 5 2 5 5 95*x *y - 53*x *y - 78*x *y + 69*x *y + 58*x *y - 58*x 4 5 4 4 4 3 4 2 4 + 64*x *y + 93*x *y - 21*x *y + 24*x *y - 13*x *y 4 3 5 3 4 3 3 3 2 - 28*x - 57*x *y - 78*x *y - 44*x *y + 37*x *y 3 3 2 5 2 4 2 3 2 2 - 64*x *y - 95*x - 71*x *y - 69*x *y - x *y - 49*x *y 2 2 5 4 3 2 + 77*x *y + 48*x + 38*x*y + 93*x*y - 65*x*y - 83*x*y 5 4 3 2 + 25*x*y + 51*x + 35*y - 18*y - 59*y + 73*y - y + 31 % A more conventional alternative is ... % procedure r; randpoly(y)$ randpoly({x}, coeffs = r); % or, in fact, equivalently ... % randpoly({x}, coeffs = procedure r; randpoly(y)); randpoly({x, y}, dense); 5 4 4 3 2 3 3 85*x + 43*x *y + 68*x + 87*x *y - 93*x *y - 20*x 2 2 2 2 4 3 2 - 74*x *y - 29*x *y + 7*x + 10*x*y + 62*x*y - 86*x*y 5 4 3 2 + 15*x*y - 97*x - 53*y + 71*y - 46*y - 28*y + 79*y + 44 \end{verbatim} \appendix \newfont{\SYM}{msbm10 scaled\magstephalf} % AMS "blackboard bold" etc \newcommand{\N}{\mbox{\SYM N}} %%% {{\bf N}} \newcommand{\th}{\mbox{$^{\it th}$}} \newtheorem{prop}{Proposition} \newenvironment{proof}% {\par\addvspace\baselineskip\noindent{\bf Proof~}}% {\hspace*{\fill}$\Box$\par\addvspace\baselineskip} \section{Algorithmic background} The only part of this package that involves any mathematics that is not completely trivial is the procedure to generate a sparse set of monomials of specified maximum and minimum total degrees in a specified set of variables. This involves some combinatorics, and the Maple implementation calls some procedures from the Maple Combinatorial Functions Package {\tt combinat} (of which I have implemented restricted versions in REDUCE). Given the maximum possible number $N$ of terms (in a dense polynomial), the required number of terms (in the sparse polynomial) is selected as a random subset of the natural numbers up to $N$, where each number indexes a term. In the univariate case these indices are used directly as monomial exponents, but in the multivariate case they are converted to monomial exponent vectors using a lexicographic ordering. \subsection{Numbers of polynomial terms} By explicitly enumerating cases with 1, 2, etc.\ variables, as indicated by the inductive proof below, one deduces that: \begin{prop} In $n$ variables, the number of distinct monomials having total degree precisely $r$ is $^{r+n-1}C_{n-1}$, and the maximum number of distinct monomials in a polynomial of maximum total degree $d$ is $^{d+n}C_n$. \end{prop} \begin{proof} Suppose the first part of the proposition is true, namely that there are at most \[ N_h(n,r) = {}^{r+n-1}C_{n-1} \] distinct monomials in an $n$-variable {\em homogeneous\/} polynomial of total degree $r$. Then there are at most \[ N(d,r) = \sum_{r=0}^d {}^{r+n-1}C_{n-1} = {}^{d+n}C_n \] distinct monomials in an $n$-variable polynomial of maximum total degree $d$. The sum follows from the fact that \[ {}^{r+n}C_n = \frac{(r+n)^{\underline n}}{n!} \] where $x^{\underline n} = x(x-1)(x-2)\cdots(x-n+1)$ denotes a falling factorial, and \[ \sum_{a \leq x < b} x^{\underline n} = \left. \frac{x^{\underline{n+1}}}{n+1} \right|_a^b. \] (See, for example, D. H. Greene \& D. E. Knuth, {\it Mathematics for the Analysis of Algorithms}, Birkh\"auser, Second Edn.\ 1982, equation (1.37)). Hence the second part of the proposition follows from the first. The proposition holds for 1 variable ($n = 1$), because there is clearly 1 distinct monomial of each degree precisely $r$ and hence at most $d+1$ distinct monomials in a polynomial of maximum degree $d$. Suppose that the proposition holds for $n$ variables, which are represented by the vector $X$. Then a homogeneous polynomial of degree $r$ in the $n+1$ variables $X$ together with the single variable $x$ has the form \[ x^r P_0(X) + x^{r-1} P_1(X) + \cdots + x^0 P_r(X) \] where $P_s(X)$ represents a polynomial of maximum total degree $s$ in the $n$ variables $X$, which therefore contains at most $^{s+n}C_n$ distinct monomials. The homogeneous polynomial of degree $r$ in $n+1$ terms therefore contains at most \[ \sum_{s=0}^r {}^{s+n}C_n = {}^{r+n+1}C_{n+1} \] distinct monomials. Hence the proposition holds for $n+1$ variables, and therefore by induction it holds for all $n$. \end{proof} \subsection{Mapping indices to exponent vectors} The previous proposition is also the basis of the algorithm to map term indices $m \in \N$ to exponent vectors $v \in \N^n$, where $n$ is the number of variables. Define a norm $\|\cdot\|$ on exponent vectors by $\|v\| = \sum_{i=1}^n v_i$, which corresponds to the total degree of the monomial. Then, from the previous proposition, the number of exponent vectors of length $n$ with norm $\|v\| \leq d$ is $N(n,d) = {}^{d+n}C_n$. The elements of the $m\th$ exponent vector are constructed recursively by applying the algorithm to successive tail vectors, so let a subscript denote the length of the vector to which a symbol refers. The aim is to compute the vector of length $n$ with index $m = m_n$. If this vector has norm $d_n$ then the index and norm must satisfy \[ N(n,d_n-1) \leq m_n < N(n,d_n), \] which can be used (as explained below) to compute $d_n$ given $n$ and $m_n$. Since there are $N(n,d_n-1)$ vectors with norm less than $d_n$, the index of the $(n-1)$-element tail vector must be given by $m_{n-1} = m_n - N(n,d_n-1)$, which can be used recursively to compute the norm $d_{n-1}$ of the tail vector. From this, the first element of the exponent vector is given by $v_1 = d_n - d_{n-1}$. The algorithm therefore has a natural recursive structure that computes the norm of each tail subvector as the recursion stack is built up, but can only compute the first term of each tail subvector as the recursion stack is unwound. Hence, it constructs the exponent vector from right to left, whilst being applied to the elements from left to right. The recursion is terminated by the observation that $v_1 = d_1 = m_1$ for an exponent vector of length $n = 1$. The main sub-procedure, given the required length $n$ and index $m_n$ of an exponent vector, must return its norm $d_n$ and the index of its tail subvector of length $n-1$. Within this procedure, $N(n,d)$ can be efficiently computed for values of $d$ increasing from 0, for which $N(n,0) = {}^nC_n = 1$, until $N(n,d) > m$ by using the observation that \[ N(n,d) = {}^{d+n}C_n = \frac{(d+n)(d-1+n)\cdots(1+n)}{d!}. \] \end{document}