File: algorithms-multivar.chapt.txt

package info (click to toggle)
yacas 1.3.6-2
  • links: PTS
  • area: main
  • in suites: bullseye, buster, sid, stretch
  • size: 7,176 kB
  • ctags: 3,520
  • sloc: cpp: 13,960; java: 12,602; sh: 11,401; makefile: 552; perl: 517; ansic: 381
file content (370 lines) | stat: -rw-r--r-- 15,420 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
		Sparse representations

	    The sparse tree data structure

Yacas has a sparse tree object for use as a storage for storing 
(key,value) pairs for which the following properties hold:

*	(key, value1) + (key, value2) = (key, value1+value2)

In addition, for multiplication the following rule is obeyed:

*	(key1, value1) * (key2, value2) = (key1+key2, value1*value2)

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.

	    Internal structure

A key is defined to be a list of integer numbers ($ n[1] $, ..., $ n[m] $).
Thus for a two-dimensional key, one item in the sparse tree database
could be reflected as the (key,value) pair { {{1,2},3} }, which states that 
element {(1,2)} has value {3}. (Note: this is not the way it is stored
in the database!).

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}}}};
{CreateSparseTree} makes a database with exactly one item.
Items can now be obtained from the sparse tree with {SparseTreeGet}.

	In> SparseTreeGet({1,2},r)
	Out> 3;
	In> SparseTreeGet({1,3},r)
	Out> 0;
And values can also be set or changed:

	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 variable {Current} represents the current value, and can
be used to determine the new value. {SparseTreeSet} destructively
modifies the original, and returns the new value. If the key
pair was not found, it is added to the tree.

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[i][j] $ 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[i][j] $ are defined by a function
$ f(i,j) $ (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:
*	{CreateSparseTree(coefs,fact)} - 
Create a sparse tree with one monomial, where 'coefs' is the 
key, and 'fact' the value. 'coefs' should be a list of integers.
*	{SparseTreeMap(op,depth,tree)} - 
Walk over the sparse tree, one element at a time, and apply the 
function "op" on the arguments (key,value). The 'value' in the
tree is replaced by the value returned by the {op} function. 'depth'
signifies the dimension of the tree (number of indices in the key).
*	{SparseTreeScan(op,depth,tree)} - 
Same as SparseTreeMap, but without changing elements.
*	{AddSparseTrees(depth,x,y)}, {MultiplyAddSparseTrees(depth,x,y,coefs,fact)} -
Add sparse tree 'y' to sparse tree 'x', destructively.
in the {MultiplyAdd} case, the monomials are treated as 
if they were multiplied by a monomial with coefficients 
with the (key,value) pair (coefs,fact). 'depth'
signifies the dimension of the tree (number of indices in the key).
*	{SparseTreeGet(key,tree)} - 
return value stored for key in the tree.
*	{SparseTreeSet(key,tree,newvalue)} - 
change the value stored for the key to newvalue. If the key was not found then
{newvalue} is stored as a new item. The variable {Current} is set to the
old value (or zero if the key didn't exist in the tree) before evaluating
{newvalue}. 


		Implementation of multivariate polynomials

This section describes the implementation of multivariate
polynomials in Yacas.

Concepts and ideas are taken from the books [Davenport <i>et al.</i> 1989] and [von zur Gathen <i>et al.</i> 1999].


	    Definitions

The following definitions define multivariate polynomials,
and the functions defined on them that are of interest
for using such multivariates.

A <i>term</i> is an object which can be written as 

$$ c*x[1]^n[1]*x[2]^n[2]* ... *x[m]^n[m] $$
for $ m $ variables ($ x[1] $, ..., $ x[m] $). The numbers $ n[m] $
are integers. $ c $ is called a <i>coefficient</i>, and 
$ x[1]^n[1]*x[2]^n[2]* ... *x[m]^n[m] $ a <i>monomial</i>.

A <i>multivariate polynomial</i> 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 <i>coefficient</i> 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 <i>lexicographic order</i> has 
been chosen:

*	first an ordering of variables is chosen, ( $ x[1] $, ..., $ x[m] $)
*	for the exponents of a monomial, $ a $ = ($ a[1] $, ..., $ a[m] $) the
lexicographic order first looks at the first exponent, $ a[1] $,
to determine which of the two monomials comes first in the multivariate.
If the two exponents are the same, the next exponent is considered.

This method is called <i>lexicographic</i> 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 <i>multivariate polynomial</i>

$$ f = Sum(a,a[max],a[min],c[a]*x^a) $$

with a monomial order:

*	1. $ c[a]*x^a $ is a <i>term</i> of the multivariate.
*	1. the <i>multidegree</i> of $ f $ is $ mdeg(f) := a[max] $.
*	1. the <i>leading coefficient</i> of $ f $ is $ lc(f):=c[mdeg(f)] $,
for the first term with non-zero coefficient.
*	1. the <i>leading monomial</i> of $ f $ is $ lm(f):=x^mdeg(f) $.
*	1. the <i>leading term</i> of $ f $ is $ lt(f):=lc(f)*lm(f) $.

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.
$ mdeg(f) $ can be thought of as a 'key', and $ lc(f) $ as a 'value'.


The <i>multicontent</i>, $ multicont(f) $, is defined to be a term that
divides all the terms in $ f $, and is the term described by
($ Min(a) $, $ Gcd(c) $), with $ Gcd(c) $ the GCD of all the coefficients,
and $ Min(a) $ the lowest exponents for each variable, occurring
in $ f $ for which $ c $ is non-zero.

The <i>multiprimitive part</i> is then defined as $ pp(f):=f$/$multicont(f) $.

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:
$$ muadd(f,t,g) := f+t*g $$.
This allows for both adding two polynomials ($ t:=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.


	    Representation

For the representation of polynomials, on computers it is natural to do this
in an array: ($a[1]$, $a[2]$, ..., $a[n]$) for a univariate polynomial, and the equivalent
for multivariates. This is called a <i>dense</i> representation, because all
the coefficients are stored, even if they are zero. 
Computers are efficient at dealing with arrays. However, in the case of 
multivariate polynomials, arrays can become
rather large, requiring a lot of storage and processing power even
to add two such polynomials. For instance, $ x^200*y^100*z^300+1 $
could take 6000000 places in an array for the coefficients. Of course
variables could be substituted for the single factors, $ p:=x^200 $
etc., but it requires an additional ad hoc step.

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}}}});
The first item in the main list is {{2,{{0,3}}}}, which states that
there is a term of the form $ x^2*y^0*3 $. The second item 
states that there are two terms, $ x^0*y^1*1 $ and $ x^0*y^0*0 = 0 $.

This representation is sparse:

	In> r:=MM(x^1000+x)
	Out> MultiNomial({x},{{1000,1},{1,1}});
and allows for easy multiplication:

	In> r*r
	Out> MultiNomial({x},{{2000,1},{1001,2},
	  {2,1},{0,0}});
	In> NormalForm(%)
	Out> x^2000+2*x^1001+x^2;


	    Internal code organization

The implementation of multivariates can be divided in three levels.

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");
where "multivar.rep/sparsenomial.ys" is the file implementing the 
driver (this is also the default driver, so the above command would
not change any thing).

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.

	    Driver interface

The driver should implement the following interface:

*	{CreateTerm(vars,{exp,coef})} - 
create a multivariate polynomial with one term, in the
variables defined in 'var', with the (key,value) pair
(coefs,fact)
*	{MultiNomialAdd(multi1, multi2) } -
add two multivars, and (possibly) destructively modify multi1 to contain
the result: [ multi1 := multi1 + multi2; multi1; ];
*	{MultiNomialMultiplyAdd(multi1, multi2,exp,coef)} - 
add two multivars, and (possibly) destructively modify multi1 to contain
the result. multi2 is considered multiplied by a term
represented by the (key,value) pair (exp,coef).
[ multi1 := multi1 + term * multi2; multi1; ];
*	{MultiNomialNegate(multi)} - 
negate a multivar, returning -multi, and destructively changing
the original. [ multi := - multi; multi1; ];
*	{MultiNomialMultiply(multi1,multi2)} - 
Multiply two multivars, and (possibly) destructively modify multi1 to contain
the result, returning the result: [ multi1 := multi1 * multi2; multi1; ];
*	{NormalForm(multi)} - 
convert MultiNomial to normal form (as would be typed in be the user).
This is part of the driver because the driver might be able to do
this more efficiently than code above it which can use ScanMultiNomial.
*	{MultiLeadingTerm(multi)} - 
return the (key,value) pair (mdeg(f),lc(f)) representing the leading
term. This is all the information needed about the leading term, 
and thus the leading coefficient and multidegree can be extracted from it.
*	{MultiDropLeadingZeroes(multi)} - 
remove leading terms with zero factors.
*	{MultiTermLess(x,y)} - 
for two (key,value) pairs, return {True} if $x<y$, where the operation {<} is the one used for
the representation, and {False} otherwise.
*	{ScanMultiNomial(op,multi)} - 
traverse all the terms of the multivariate, applying the
function 'op' to each (key,value) pair (exp,coef). The monomials
are traversed in the ordering defined by MultiTermLess. 'op' should
be a function accepting two arguments.
*	{MultiZero(multi)} - 
return {True} if the multivariate is zero (all coefficients are zero),
{False} otherwise.