In addition, for multiplication the following rule is obeyed:
The last is optional. For multivariate polynomials (described elsewhere) both hold, but for matrices, only the addition property holds. The function MultiplyAddSparseTrees (described below) should not be used in these cases.
The storage is recursive. The sparse tree begins with a list of objects {n1,tree1} for values of n1 for the first item in the key. The tree1 part then contains a sub-tree for all the items in the database for which the value of the first item in the key is n1.
The above single element could be created with
In> r:=CreateSparseTree({1,2},3)
Out> {{1,{{2,3}}}};
|
In> SparseTreeGet({1,2},r)
Out> 3;
In> SparseTreeGet({1,3},r)
Out> 0;
|
In> SparseTreeSet({1,2},r,Current+5)
Out> 8;
In> r
Out> {{1,{{2,8}}}};
In> SparseTreeSet({1,3},r,Current+5)
Out> 5;
In> r
Out> {{1,{{3,5},{2,8}}}};
|
The sparse tree can be traversed, one element at a time, with SparseTreeScan:
In> SparseTreeScan(Hold({{k,v},Echo({k,v})}),2,r)
{1,3} 5
{1,2} 8
|
An example of the use of this function could be multiplying a sparse matrix with a sparse vector, where the entire matrix can be scanned with SparseTreeScan, and each non-zero matrix element $A _{\left( i, j\right) }$ can then be multiplied with a vector element $v _{j}$, and the result added to a sparse vector $w _{i}$, using the SparseTreeGet and SparseTreeSet functions. Multiplying two sparse matrices would require two nested calls to SparseTreeScan to multiply every item from one matrix with an element from the other, and add it to the appropriate element in the resulting sparse matrix.
When the matrix elements $A _{\left( i, j\right) }$ are defined by a function $f\left( i, j\right) $ (which can be considered a dense representation), and it needs to be multiplied with a sparse vector $v _{j}$, it is better to iterate over the sparse vector $v _{j}$. Representation defines the most efficient algorithm to use in this case.
The API to sparse trees is:
Concepts and ideas are taken from the books [Davenport et al. 1989] and [von zur Gathen et al. 1999].
A term is an object which can be written as
$$c x _{1} ^{n _{1}} x _{2} ^{n _{2}} \mathrm{ ... } x _{m} ^{n _{m}}$$ for $m$ variables ($x _{1}$, ..., $x _{m}$). The numbers $n _{m}$ are integers. $c$ is called a coefficient, and $x _{1} ^{n _{1}} x _{2} ^{n _{2}} \mathrm{ ... } x _{m} ^{n _{m}}$ a monomial.
A multivariate polynomial is taken to be a sum over terms.
We write $c _{a} x ^{a}$ for a term, where $a$ is a list of powers for the monomial, and $c _{a}$ the coefficient of the term.
It is useful to define an ordering of monomials, to be able to determine a canonical form of a multivariate.
For the currently implemented code the lexicographic order has been chosen:
This method is called lexicographic because it is similar to the way words are ordered in a usual dictionary.
For all algorithms (including division) there is some freedom in the ordering of monomials. One interesting advantage of the lexicographic order is that it can be implemented with a recursive data structure, where the first variable, $x _{1}$ can be treated as the main variable, thus presenting it as a univariate polynomial in $x _{1}$ with all its terms grouped together.
Other orderings can be used, by re-implementing a part of the code dealing with multivariate polynomials, and then selecting the new code to be used as a driver, as will be described later on.
Given the above ordering, the following definitions can be stated:
For a non-zero multivariate polynomial
$$f = \sum _{a = a _{\mathrm{ max }}} ^{a _{\mathrm{ min }}} c _{a} x ^{a}$$
with a monomial order:
The above define access to the leading monomial, which is used for divisions, gcd calculations and the like. Thus an implementation needs be able to determine {mdeg(f),lc(f)} . Note the similarity with the (key,value) pairs described in the sparse tree section. $\mathrm{ mdeg }\left( f\right) $ can be thought of as a 'key', and $lc\left( f\right) $ as a 'value'.
The multicontent, $\mathrm{ multicont }\left( f\right) $, is defined to be a term that divides all the terms in $f$, and is the term described by ($\min \left( a\right) $, $\mathrm{ Gcd }\left( c\right) $), with $\mathrm{ Gcd }\left( c\right) $ the GCD of all the coefficients, and $\min \left( a\right) $ the lowest exponents for each variable, occurring in $f$ for which $c$ is non-zero.
The multiprimitive part is then defined as $pp\left( f\right) \equiv f$/$\mathrm{ multicont }\left( f\right) $.
For a multivariate polynomial, the obvious addition and (distributive) multiplication rules hold:
(a+b) + (c+d) := a+b+c+d |
a*(b+c) := (a*b)+(a*c) |
These are supported in the Yacas system through a multiply-add operation: $$\mathrm{ muadd }\left( f, t, g\right) \equiv f + t g.$$ This allows for both adding two polynomials ($t\equiv 1$), or multiplication of two polynomials by scanning one polynomial, and multiplying each term of the scanned polynomial with the other polynomial, and adding the result to the polynomial that will be returned. Thus there should be an efficient muadd operation in the system.
An alternative is to store only the terms for which the coefficients are non-zero. This adds a little overhead to polynomials that could efficiently be stored in a dense representation, but it is still little memory, whereas large sparse polynomials are stored in acceptable memory too. It is of importance to still be able to add, multiply divide and get the leading term of a multivariate polynomial, when the polynomial is stored in a sparse representation.
For the representation, the data structure containing the (exponents,coefficient) pair can be viewed as a database holding (key,value) pairs, where the list of exponents is the key, and the coefficient of the term is the value stored for that key. Thus, for a variable set {x,y} the list {{1,2},3} represents $3 x y ^{2}$.
Yacas stores multivariates internally as MultiNomial (vars, terms), where vars is the ordered list of variables, and terms some object storing all the (key, value) pairs representing the terms. Note we keep the storage vague: the terms placeholder is implemented by other code, as a database of terms. The specific representation can be configured at startup (this is described in more detail below).
For the current version, Yacas uses the 'sparse tree' representation, which is a recursive sparse representation. For example, for a variable set {x,y,z}, the 'terms' object contains a list of objects of form {deg,terms}, one for each degree deg for the variable 'x' occurring in the polynomial. The 'terms' part of this object is then a sub-sparse tree for the variables {y,z}.
An explicit example:
In> MM(3*x^2+y)
Out> MultiNomial({x,y},{{2,{{0,3}}},{0,{{1,1},
{0,0}}}});
|
This representation is sparse:
In> r:=MM(x^1000+x)
Out> MultiNomial({x},{{1000,1},{1,1}});
|
In> r*r
Out> MultiNomial({x},{{2000,1},{1001,2},
{2,1},{0,0}});
In> NormalForm(%)
Out> x^2000+2*x^1001+x^2;
|
At the top level are the routines callable by the user or the rest of the system: MultiDegree, MultiDivide, MultiGcd, Groebner, etc. In general, this is the level implementing the operations actually desired.
The middle level does the book-keeping of the MultiNomial(vars,terms) expressions, using the functionality offered by the lowest level.
For the current system, the middle level is in multivar.rep/ sparsenomial.ys, and it uses the sparse tree representation implemented in sparsetree.ys.
The middle level is called the 'driver', and can be changed, or re-implemented if necessary. For instance, in case calculations need to be done for which dense representations are actually acceptable, one could write C++ implementing above-mentioned database structure, and then write a middle-level driver using the code. The driver can then be selected at startup. In the file 'yacasinit.ys' the default driver is chosen, but this can be overridden in the .yacasrc file or some file that is loaded, or at the command line, as long as it is done before the multivariates module is loaded (which loads the selected driver). Driver selection is as simple as setting a global variable to contain a file name of the file implementing the driver:
Set(MultiNomialDriver, "multivar.rep/sparsenomial.ys"); |
The choice was made for static configuration of the driver before the system starts up because it is expected that there will in general be one best way of doing it, given a certain system with a certain set of libraries installed on the operating system, and for a specific version of Yacas. The best version can then be selected at start up, as a configuration step. The advantage of static selection is that no overhead is imposed: there is no performance penalty for the abstraction layers between the three levels.
Integrate can have its own set of rules for specific integrals, which might return a correct answer immediately. Alternatively, it calls the function AntiDeriv, to see if the anti-derivative can be determined for the integral requested. If this is the case, the anti-derivative is used to compose the output.
If the integration algorithm cannot perform the integral, the expression is returned unsimplified.
For the purpose of setting up the integration table, a few declaration functions have been defined, which use some generalized pattern matchers to be more flexible in recognizing expressions that are integrable.
The calling sequence for IntFunc is
IntFunc(variable,pattern,antiderivative) |
For instance, for the function $\cos x$ there is a declaration:
IntFunc(x,Cos(_x),Sin(x)); |
The fact that the second argument is a pattern means that each occurrence of the variable to be matched should be referred to as _x, as in the example above.
IntFunc generalizes the integration implicitly, in that it will set up the system to actually recognize expressions of the form $\cos \left( a x + b\right) $, and return $\frac{\sin \left( a x + b\right) }{a} $ automatically. This means that the variables $a$ and $b$ are reserved, and can not be used in the pattern. Also, the variable used (in this case, _x is actually matched to the expression passed in to the function, and the variable var is the real variable being integrated over. To clarify: suppose the user wants to integrate $\cos \left( c y + d\right) $ over $y$, then the following variables are set:
When functions are multiplied by constants, that situation is handled by the integration rule that can deal with univariate polynomials multiplied by functions, as a constant is a polynomial of degree zero.
The general form for declaring anti-derivatives for such expressions is:
IntPureSquare(variable, pattern, sign2, sign0, antiderivative) |
The expression is searched for the pattern, where the variable can match to a sub-expression of the form $a x ^{2} + b$, and for which both $a$ and $b$ are numbers and $a \mathrm{ sign2 } > 0$ and $b \mathrm{ sign0 } > 0$.
As an example:
IntPureSquare(x,num_IsFreeOf(var)/(_x), 1,1,(num/(a*Sqrt(b/a)))* ArcTan(var/Sqrt(b/a))); |
It has been setup much like the integration algorithm. If the transformation algorithm cannot perform the transform, the expression (in theory) is returned unsimplified. Some cases may still erroneously return Undefined or Infinity.
The last operational property dealing with integration is not yet fully bug-tested, it sometimes returns Undefined or Infinity if the integral returns such.
For the purpose of setting up the transform table, a few declaration functions have been defined, which use some generalized pattern matchers to be more flexible in recognizing expressions that are transformable.
The calling sequence for LapTranDef is
LapTranDef( in, out ) |
Currently in must be a variable of _t and out must be a function of s. For instance, for the function $\cos t$ there is a declaration:
LapTranDef( Cos(_t), s/(s^2+1) ); |
The fact that the first argument is a pattern means that each occurrence of the variable to be matched should be referred to as _t, as in the example above.
LapTranDef generalizes the transform implicitly, in that it will set up the system to actually recognize expressions of the form $\cos a t$ and $\cos \frac{t}{a} $ , and return the appropriate answer. The way this is done is by three separate rules for case of t itself, a*t and t/a. This is similar to the MatchLinear function that Integrate uses, except LaplaceTransforms must have b=0.
Without loss of generality, the coefficients $a _{i}$ of a polynomial $$p = a _{n} x ^{n} + \mathrm{ ... } + a _{0}$$ can be considered to be rational numbers, as real-valued numbers are truncated in practice, when doing calculations on a computer.
Assuming that the leading coefficient $a _{n} = 1$, the polynomial $p$ can also be written as $$p = p _{1} ^{n _{1}} \mathrm{ ... } p _{m} ^{n _{m}},$$ where $p _{i}$ are the $m$ distinct irreducible monic factors of the form $p _{i} = x - x _{i}$, and $n _{i}$ are multiplicities of the factors. Here the roots are $x _{i}$ and some of them may be complex. However, complex roots of a polynomial with real coefficients always come in conjugate pairs, so the corresponding irreducible factors should be taken as $p _{i} = x ^{2} + c _{i} x + d _{i}$. In this case, there will be less than $m$ irreducible factors, and all coefficients will be real.
To find roots, it is useful to first remove the multiplicities, i.e. to convert the polynomial to one with multiplicity 1 for all irreducible factors, i.e. find the polynomial $p _{1} \mathrm{ ... } p _{m}$. This is called the "square-free part" of the original polynomial $p$.
The square-free part of the polynomial $p$ can be easily found using the polynomial GCD algorithm. The derivative of a polynomial $p$ can be written as: $$p' = \sum _{i = 1} ^{m} p _{1} ^{n _{1}} \mathrm{ ... } n _{i} p _{i} ^{n _{i} - 1} \left( \frac{\partial}{\partial x}p _{i}\right) \mathrm{ ... } p _{m} ^{n _{m}}.$$
The g.c.d. of $p$ and $p'$ equals $$\mathrm{ Gcd }\left( p, p'\right) = \prod _{i = 1} ^{m} p _{i} ^{n _{i} - 1}.$$ So if we divide $p$ by $\mathrm{ Gcd }\left( p, p'\right) $, we get the square-free part of the polynomial: $$\mathrm{ SquareFree }\left( p\right) \equiv \mathrm{ Div }\left( p, \mathrm{ Gcd }\left( p, p'\right) \right) = p _{1} \mathrm{ ... } p _{m}.$$
In what follows we shall assume that all polynomials are square-free with rational coefficients. Given any polynomial, we can apply the functions SquareFree and Rationalize and reduce it to this form. The function Rationalize converts all numbers in an expression to rational numbers. The function SquareFree returns the square-free part of a polynomial. For example:
In> Expand((x+1.5)^5) Out> x^5+7.5*x^4+22.5*x^3+33.75*x^2+25.3125*x +7.59375; In> SquareFree(Rationalize(%)) Out> x/5+3/10; In> Simplify(%*5) Out> (2*x+3)/2; In> Expand(%) Out> x+3/2; |
$$p _{0} = p\left( x\right) ,$$ $$p _{1} = p'\left( x\right) ,$$ $$p _{i} = - \mathrm{ remainder }\left( p _{i - 2}, p _{i - 1}\right) ,$$ where $\mathrm{ remainder }\left( p, q\right) $ is the remainder of division of polynomials $p$, $q$.
The polynomial $p$ can be assumed to have no multiple factors, and thus $p$ and $p'$ are relatively prime. The sequence of polynomials in the Sturm sequence are (up to a minus sign) the consecutive polynomials generated by Euclid's algorithm for the calculation of a greatest common divisor for $p$ and $p'$, so the last polynomial $p _{n}$ will be a constant.
In Yacas, the function SturmSequence(p) returns the Sturm sequence of $p$, assuming $p$ is a univariate polynomial in $x$, $p = p\left( x\right) $.
Given a Sturm sequence $S = \mathrm{ SturmSequence }\left( p\right) $ of a polynomial $p$, the variation in the Sturm sequence $V\left( S, y\right) $ is the number of sign changes in the sequence $p _{0}$, $p _{1}$ , ... , $p _{n}$, evaluated at point $y$, and disregarding zeroes in the sequence.
Sturm's theorem states that if $a$ and $b$ are two real numbers which are not roots of $p$, and $a < b$, then the number of roots between $a$ and $b$ is $V\left( S, a\right) - V\left( S, b\right) $. A proof can be found in Knuth, The Art of Computer Programming, Volume 2, Seminumerical Algorithms.
For $a$ and $b$, the values $ - \infty $ and $\infty $ can also be used. In these cases, $V\left( S, \infty \right) $ is the number of sign changes between the leading coefficients of the elements of the Sturm sequence, and $V\left( S, - \infty \right) $ the same, but with a minus sign for the leading coefficients for which the degree is odd.
Thus, the number of real roots of a polynomial is $V\left( S, - \infty \right) - V\left( S, \infty \right) $. The function NumRealRoots(p) returns the number of real roots of $p$.
The function SturmVariations(S,y) returns the number of sign changes between the elements in the Sturm sequence $S$, at point $x = y$:
In> p:=x^2-1
Out> x^2-1;
In> S:=SturmSequence(p)
Out> {x^2-1,2*x,1};
In> SturmVariations(S,-Infinity)- \
SturmVariations(S,Infinity)
Out> 2;
In> NumRealRoots(p)
Out> 2;
In> p:=x^2+1
Out> x^2+1;
In> S:=SturmSequence(p)
Out> {x^2+1,2*x,-1};
In> SturmVariations(S,-Infinity)- \
SturmVariations(S,Infinity)
Out> 0;
In> NumRealRoots(p)
Out> 0;
|
$$\left| a\right| \leq \max \left( \left| \frac{a _{n - 1}}{a _{n}} \right| , \sqrt{\left| \frac{a _{n - 2}}{a _{n}} \right| }, \mathrm{ ... }, \sqrt[n]{\left| \frac{a _{0}}{a _{n}} \right| }\right) ,$$ for all real roots $a$ of $p$. This gives the upper bound on the absolute value of the roots of the polynomial in question. if $p\left( 0\right) \neq 0$, the minimum bound can be obtained also by considering the upper bound of $p\left( \frac{1}{x} \right) x ^{n}$, and taking $\frac{1}{\mathrm{ bound }} $.
We thus know that given $$a _{\mathrm{ max }} = \mathrm{ MaximumBound }\left( p\right) $$ and $$a _{\mathrm{ min }} = \mathrm{ MinimumBound }\left( p\right) $$ for all roots $a$ of polynomial, either $$ - a _{\mathrm{ max }}\leq a\leq - a _{\mathrm{ min }}$$ or $$a _{\mathrm{ min }}\leq a\leq a _{\mathrm{ max }}.$$
Now we can start the search for the bounds on all roots. The search starts with initial upper and lower bounds on ranges, subdividing ranges until a range contains only one root, and adding that range to the resulting list of bounds. If, when dividing a range, the middle of the range lands on a root, care must be taken, because the bounds should not be on a root themselves. This can be solved by observing that if $c$ is a root, $p$ contains a factor $x - c$, and thus taking $p\left( x + c\right) $ results in a polynomial with all the roots shifted by a constant $ - c$, and the root $c$ moved to zero, e.g. $p\left( x + c\right) $ contains a factor $x$. Thus a new ranges to the left and right of $c$ can be determined by first calculating the minimum bound $M$ of $\frac{p\left( x + c\right) }{x} $. When the original range was ($a$,$b$), and $c = \frac{a + b}{2} $ is a root, the new ranges should become ($a$,$c - M$) and ($c + M$,$b$).
In Yacas, MimimumBound(p) returns the lower bound described above, and MaximumBound(p) returns the upper bound on the roots in $p$. These bounds are returned as rational numbers. BoundRealRoots(p) returns a list with sublists with the bounds on the roots of a polynomial:
In> p:=(x+20)*(x+10)
Out> (x+20)*(x+10);
In> MinimumBound(p)
Out> 10/3;
In> MaximumBound(p)
Out> 60;
In> BoundRealRoots(p)
Out> {{-95/3,-35/2},{-35/2,-10/3}};
In> N(%)
Out> {{-31.6666666666,-17.5},
{-17.5,-3.3333333333}};
|
It should be noted that since all calculations are done with rational numbers, the algorithm for bounding the roots is very robust. This is important, as the roots can be very unstable for small variations in the coefficients of the polynomial in question (see Davenport).
The bisection method is more robust, but slower. It works by taking the middle of the range, and checking signs of the polynomial to select the half-range where the root is. As there is only one root in the range ($a$,$b$), in general it will be true that $p\left( a\right) p\left( b\right) < 0$, which is assumed by this method.
Yacas finds the roots by first trying the secant method, starting in the middle of the range, $c = \frac{a + b}{2} $. If this fails the bisection method is tried.
The function call to find the real roots of a polynomial $p$ in variable $x$ is FindRealRoots(p), for example:
In> p:=Expand((x+3.1)*(x-6.23))
Out> x^2-3.13*x-19.313;
In> FindRealRoots(p)
Out> {-3.1,6.23};
In> p:=Expand((x+3.1)^3*(x-6.23))
Out> x^4+3.07*x^3-29.109*x^2-149.8199\
In> *x-185.59793;
In> p:=SquareFree(Rationalize( \
In> Expand((x+3.1)^3*(x-6.23))))
Out> (-160000*x^2+500800*x+3090080)/2611467;
In> FindRealRoots(p)
Out> {-3.1,6.23};
|