POLY : A new polynomial data
structure for Maple 17 ?
Michael Monagan and Roman Pearce

Abstract
We demonstrate how a new data structure for sparse distributed polynomials in the Maple 17 kernel significantly accelerates several key Maple library routines. The POLY data structure and its associated kernel operations
are programmed for compactness, scalability, and low overhead. This allows
polynomials to have tens of millions of terms, increases parallel speedup, andimproves the performance of Maple library routines.

1 Introduction
Figure 1 shows the default polynomial data structure in Maple 16 and previous versions of Maple for the polynomial f = 9xy 3 z −4y 3 z 2 −6xy 2 z −8x3 −5.
It is a “sum-of-products” where each monomial is stored a separate Maple
object, a PROD. To compute the degree of f , a coefficient in x, test for asubexpression, or do almost anything else, the Maple kernel must recursively
descend through multiple levels of dags. This involves extensive branching
and random memory access, which prevents Maple from achieving highperformance on modern computer processors.
For example, to compute the degree of f in x, when Maple sees x in
the monomial PROD 7 x 1 y 3 z 1 , because Maple does not have a dedicatedpolynomial data structure, Maple does not know that the next factor does
not have x in it. The monomial could be x(x + y)3 z or x sin(x)3 z. Thus,
Department of Mathematics, Simon Fraser University, Burnaby, British Columbia, V5A
1S6, Canada. e-mail: mmonagan@cecm.sfu.ca and e-mail: rpearcea@cecm.sfu.ca? This work was supported by Maplesoft and the National Science and Engineering Research Council (NSERC) of Canada.

12Michael Monagan and Roman Pearce

Fig. 1 Maple’s sum-of-products representation encodes each object in an array of words
where the first word encodes the type and length (in words) of the object.
PROD 7

SUM 11

9

1

y

3

z

1

PROD 5y3z2

PROD 7x1

y2

PROD 3

x

3

−8

−5

1

x

−4

−6

z

1

computing the degree of a polynomial with t terms in n variables is O(nt) in
Maple.
Monomial operations are particularly slow. To multiply xy 3 z by xz Maple
forms the product PROD 11 x 1 y 3 z 1 x 1 z 1 . It then “simplifies” this product to obtain PROD 7 x 2 y 3 z 2 . Since the variables in a PROD are not sorted,
Maple uses an O(n2 ) double loop to add the exponents of like variables. Finally, because Maple stores unique copies of objects, the resulting PROD object
is hashed and inserted in an internal table if it not already present. In all,there are many function calls and many loops. We estimate that Maple takes
more than 200 clock cycles to multiply monomials in 3 variables.
For comparison, Figure 2 below shows Singular’s data structure for thesame polynomial. Singular uses a linked list of terms and dense exponent
vectors to represent monomials.

Fig. 2 Singular’s polynomial representation for f = 9xy 3 z − 4y 3 z 2 − 6xy 2 z − 8x3 − 5.
POLY
x
y
z

9
1

−4
0

−6
1

−8
3

−5
0

3

3

2

0

0

1

2

1

0

0

Singular is representative of several computer algebra systems which havea dedicated distributed polynomial representation with obvious advantages.
Monomial multiplication consists of a memory allocation and a loop of fixed
length to add exponents. This is the main reason Singular is faster than
Maple at polynomial multiplication.
The relatively slow performance of Maple for polynomials lead us to develop a high performance C library for sparse polynomial arithmetic which we
integrated into Maple in [16]. To multiply two polynomials, we first convert
the Maple inputs to a special polynomial representation, then we call our li-

POLY : A new polynomial data structure for Maple 173brary to compute the product using our parallel algorithm from [18], then weconvert the result to Maple’s sum-of-products representation. This substantially improved Maple’s performance for large polynomials. However, for very
sparse polynomials, and small polynomials, the overhead of constructing thesum-of-products representation often negated the parallel speedup achievedby our library. Furthermore, algorithms in the Maple library continued to use
the sum-of-products structure, so the cost of ancillary operations like computing the degree or coefficients in a variable could easily overshadow the costof arithmetic. These problems are addressed in this paper by introducing a
new polynomial data structure in the Maple kernel. The new data structure,
called POLY, was integrated into Maple 17 which was released in April 2013.
Our paper is organized as follows. In Section 2 we describe POLY, list
some of its obvious advantages, and detail how it is integrated into Maple. In
Section 3 we describe how we implemented various Maple kernel operations
and give benchmarks demonstrating the improvement in their performance.
In Section 4 we consider the impact of the POLY data structure and improved kernel routines on Maple library codes and measure improved parallel
speedup in polynomial arithmetic. We end with a conclusion.

2 The POLY Data Structure
Figure 3 shows our new data structure for sparse distributed polynomials.
The first word is a header word, which encodes the length and type of theobject. The second word points to the variables, which are sorted in Maple’s
canonical ordering for sets. This is followed by the monomials and coefficients,
where the monomials encode the exponents and the total degree in a single
machine word. For example, for xy 2 z 3 we encode the values (6, 1, 2, 3) as theinteger 6 · 248 + 232 + 2 · 216 + 3, using 16 bits each on a 64-bit machine. Terms
are sorted into graded lexicographical order by comparing the monomials as
unsigned integers. This gives a canonical representation for the polynomial.
Small integer coefficients −262 < x < 262 are encoded as 2x + 1, so that therightmost bit is 1. Larger integers are a pointer (with rightmost bit 0) toa GMP multiprecision integer [8]. The current implementation requires all
coefficients to be integers.

Fig. 3 The new packed representation for f = 9xy 3 z − 4y 3 z 2 − 6xy 2 z − 8x3 − 5.
SEQ 4
POLY 12

5131

x

y

z

9

5032 −4

4121 −63300 −80000 −54Michael Monagan and Roman Pearce

Five advantages of the new representation are readily apparent.
1. It is much more compact. Polynomials use two words per term insteadof 2n + 3 words, where n is the number of variables. For polynomials inthree variables this saves over a factor of four.
2. Monomial comparisons become machine word comparisons and monomial multiplication becomes machine word addition (provided there is no
overflow), and monomial division becomes subtraction with a bitwise test
for failure. This dramatically speeds up polynomial arithmetic.
3. Explicitly storing variables and sorting the terms lets us perform many
common Maple idioms without looking at all of the terms, e.g. degree(f )(total degree), indets(f ) (extract the set of variables), has(f, x), andtype(f, polynom).
4. Other operations such as degree(f, x), diff(f, x), or coeff(f, x, i) (extract the coefficient of xi ) access memory sequentially to make good useof cache. We can isolate groups of exponents using masks. This eliminates
branching and loops at the level of the exponents.
5. For large polynomials, we avoid creating many small Maple objects (thePRODs) that must be simplified by Maple’s internal simplifier and storedin Maple’s simpl table, an internal hash table of all Maple objects. They
fill the simpl table and slow down Maple’s garbage collector.
The idea of packing monomials in one or more computer words is notnew; the ALTRAN computer algebra system [11] allowed the user to packmonomials in lexicographical order to conserve memory. Bruno Buchberger[4] experimented with packed monomials using the gradex lexicogrpahical ordering when he was computing Gröbner bases in the early 1960s because of
the scarcity of main memory. A number of special purpose computer algebra systems have supported monomial packing. For example, the McCaulay2
computer algebra system [9] for computing Gröbner bases supports packing
for the lexicographical and graded reverse lexicographical monomial orderings
and the Magma computer algebra system [3] dynamically packs monomials
using bytes. In [1], Bachmann and Schönemann compared the graded packing
with packings for other monomial orderings for Gröbner basis computation.
We have integrated the POLY data structure into the Maple kernel and it
is used by default as of Maple 17. When Maple parses input from the user ora file it creates the sum of products structure which is then simplified. Whenthe simplifier detects an expanded polynomial with integer coefficients and at
least two terms it creates a POLY structure when: (i) the number variablesn < β/2 for β-bit machines, (ii) the total degree d satisfies 1 < d < 2b whereb = bβ/(n + 1)c for n > 1 and b = β − 2 for n = 1, (iii) the variables are
Maple names with regular evaluation rules, e.g. {x, y1 , π} but not infinity or
undefined. Otherwise the sum-of-products format is retained. On output, andfor other cases where it was convenient, Maple recreates a sum-of-products
structure for a POLY and caches it. Thus Maple 17 uses two representations

POLY : A new polynomial data structure for Maple 17

5for polynomials, automatically selecting between the two. All conversions
between representations are automatic and invisible to the Maple user.
Note, the monomial encoding is determined solely by the number of variables in the polynomial. We use b bits for the exponent in each variable andfor the total degree. It follows that arithmetic operations between polynomials in the same variables do not need to repack monomials.
Note, for polynomials with total degree d = 1, we do not store them in asa POLY dag because Maple’s sum-of-products representation is better in thiscase. For example f = 2x + 3y + 4z + 5 is represented as SUM 9 2 x 3 y 4 z 5 1 .
This is compact as monomials are not explicitly represented. Maple’s sparselinear solver uses this polynomial representation.
We chose the graded lex ordering rather than pure lexicographical ordering for several reasons. First, the graded ordering appears more natural for
output. Second, when multiplying a × b, if the total degree d = deg a + deg b
does not overflow, that is, d < 2b , then the entire product can be computed
without overflow and with no overflow detection. This allows us to look at
only the leading terms of polynomials and predict overflow in O(1) time.
Third, in the division algorithm, if one uses pure lexicographical order, degrees in the remainder can increase and overflow the exponents. For example,
consider the following division
x2 y 5 + y 3 ÷ x2 y + xy 5in lexicographical order with x > y. The quotient is y 4 and the remainder is−xy 9 + y 3 . If we had 3 bits per variable, the y 9 would overflow. In contrast,
when a graded ordering is used the total degree of the monomials in thedivision algorithm always decreases. In our example the leading term of thedivisor would be xy 5 and the division
x2 y 5 + y 3 ÷ xy 5 + x2 ywould result in the quotient x and remainder −x3 y + y 3 .
We encode the monomial xi y j z k as d i j k where d = i + j + k. Notice,
however, that we could store d i j and recover k = d − i − j as needed toallow for more bits. We rejected this idea because it adds too much code
complexity, both in the kernel routines and in external libraries that must
support POLY. It also makes the cost of some monomial operations O(n)instead of O(1).
We expect to pack many practical problems into 64-bit words, which isthe norm for today’s computers. For example, if a polynomial has 8 variables
then we store 9 integers for each monomial using b64/9c = 7 bits each. So the
POLY dag can accommodate polynomials in 8 variables of total degree up to
127. Table 1 below shows the number of bits per exponent for a polynomial inn variables assuming a 64-bit word. Column lex shows how many bits would
be available if lexicographical order were used and one bit (at the top) wereleft unset to test for overflow. To check for overflow in a+b we would compute6Michael Monagan and Roman Pearce(a + b) ⊕ a ⊕ b. This sets the bits that are carried into (see [21]) after which
we can apply a mask.
Table 1 The number of bits per exponent for n variables in the two orderings.
n lex grlex unused n lex grlex unused
1 63
62
2 7 9
8
0
2 31
21
1 8 7
7
1
3 21
16
0 9 7
6
4
4 15
12
4 10 659
5 12
10
4 11 554
6 10
9
1 12 54
12n lex grlex unused
13 4
4
8
14 4
4
4
15 4
4
0
16 3
3
13
17 3
3
10
18 3
3
7

n lex grlex unused
19 3
3
4
20 3
3
1
21 3
2
20
22 2
2
18
... ...
...
...
31 2
2
0

Column grlex is our packing. It shows how many bits are available for eachvariable and for the total degree. For univariate polynomials we do not store
the degree because it simply duplicates the exponent. Instead we restrict thedegree to the range of Maple immediate integers (62 bits signed x is stored as
2x + 1) to avoid handling multiprecision exponents in conversions. The tablealso shows the number of unused bits at the top of each word.
When this work was presented at the Asian Symposium on Computer
Mathematics meeting in Beijing in October 2012, Joris van der Hoven asked“Why don’t you use the extra bits for increased total degree?”. For example,
in Table 1 we see that for n = 10 variables, we allocate b = 5 bits per exponent
and 5 bits for the total degree leaving 9 unused bits. The answer we gave wasthat doing so would increase the complexity of the code significantly. We
also just did not know how much difference those extra bits would make.
Over the following year, when computing determinants of various matrices ofpolynomials, we encountered three determinants which do not fit in POLY.
In all three cases, however, they would fit if we used the unused bits for total
degree. Consider the n × n Vandermonde matrix Vn where the (i, j)’th entryof Vn is xj−1
. For examplei
1 x1 x21 x31
 1 x2 x22 x32 V4 = 
 1 x3 x23 x33 
1 x4 x24 x34

The determinant Dn of Vn is a polynomial in n variables of total degree
0 + 1 + 2 + 3 + ... + n − 1 = n(n − 1)/2. Looking at Table 1, we can see that
D10 which has total degree 45 does not fit in POLY because there are only 5
bits for the total degree. However, it would easily fit if the extra 9 bits wereused for the total degree. In fact, Dn would fit for n = 11, 12, 13 and 14 as
well. This example, and the significant number of unused bits in the range
9 ≤ n ≤ 14 has convinced us that it will be worthwhile using them for the
total degree. We are presently implementing this.

POLY : A new polynomial data structure for Maple 17

7

3 Operations in the Maple Kernel
The new representation allowed us to design new high performance algorithms
for the Maple kernel. In the old data structure most operations are O(nt)where n is the number of variables and t is the number of terms. Maple
must examine the entire sum of products structure because its contents are
unknown.
With the new data structure we can often avoid doing expensive operations
on all of the terms, or we can do them much more efficiently. Our first example
is the command diff(f, x). To compute the derivative with respect to x in
the sum of products representation, Maple searches for x in each PROD. If
x is found it copies the PROD, decreases the exponent of x, and scales thecoefficient. Our algorithm for diff(f, x) first locates x in the set of variables.
If x is not there it returns zero. Otherwise let x be variable k of n, lets = b64/(n + 1)c be the width of the exponents, sk = s(n − k) is the shift
to find the exponent of x, and b = 2s − 1 is a mask of s bits. The core ofdiff(f, x) is a loop that updates the exponent of x and the total degree (ifpresent) using one subtraction.
/* subtracted from monomial */d = 1 << sk;
if (n > 1) d += 1 << (s*n);for (i=j=2; i < LENGTH(f); i+=2) {m = f[i];/* monomial */e = (m >> sk) & b; /* exponent */if (!e) continue;/* skip the constant */r[j+0] = m - d;/* new monomial */r[j+1] = mulint(f[i+1], IMMEDIATE(e), NULL);
j += 2; }The new monomials in the derivative remain sorted in the graded ordering.
For example, consider f = 9xy 3 z − 4y 3 z 2 − 6xy 2 z − 8x3 − 5. If we differentiatewith respect to x we obtain f 0 = 9y 3 z + 0 − 6y 2 z − 24x2 + 0 and the nonzero terms remain sorted in the graded monomial ordering. Thus we candifferentiate in O(n + t) instead of O(nt).
As a second example, consider degree(f, {x, z}); for f ∈ Z[x, y, z]. We
exploit a trick from Hacker’s Delight [21]. If a and b are packed monomials(a−b)⊕a⊕b sets the bits that are borrowed to subtract a−b. The absence ofunderflow in the exponent fields indicates that two monomials divide and a−bis their quotient. To test for division in {x, z} only, we construct a mask mwith 1 bits above the x and z fields and compute ((a-b) ^ a ^ b) & m. If thisis zero we know b divides a in {x, z}, so degree(b, {x, z}) < degree(a, {x, z})
and we can skip over b without computing its degree. This optimization is
also used for lcoeff which computes the leading terms of polynomials. In
the case of degree we can stop when the total degree of the monomials is
less than or equal to the largest degree already found.
Another optimization that we use is binary search. Consider has(f, x5 )and coeff(f, x, 5). In both cases it is pointless to examine terms of total

8Michael Monagan and Roman Pearce

degree less than five, but rather than test each monomial with an additional
branch in the main loop, we use binary search to adjust the endpoint for
linear search.
Our results are summarized in Table 2 which reports the speedup for kernel
operations on a polynomial with 990,000 terms in 3 variables. Note that noneof these operations are parallelized, and the cost for evaluation is added tothe other commands if you are using Maple interactively.
Table 2 Improvements for Maple kernel operations. The polynomial f is constructed asfollows f := expand( mul(randpoly(i,degree=100,dense), i=[x,y,z])):commandf;coeff(f, x, 20)
coeffs(f, x)
frontend(g, [f ])degree(f )
degree(f, x)diff(f, x)eval(f, x = 6)
expand(2xf )has(f, x101 )indets(f )lcoeff(f, x)op(f )
subs(x = y, f )taylor(f, x, 50)
type(f, polynom)

description
Maple 16
Maple 17
speedup
top level evaluation
0.162 s0.000 s→ O(n)
evaluate the variables, interactive statements incur this costcoefficient of x20
2.140 s0.005 s420x
linear search to degree 19, binary search for univariate f
extract all coefficients in x
0.979 s0.119 s8x
reorder exponents [d, x, y, z] → [x, d, y, z], sort to collect in x
subs functions for variables
3.730 s0.000 s→ O(n)
there’s nothing to do since POLY has no functions
total degree in all variables
0.175 s0.000 s→ O(1)
degree in x
0.073 s0.003 s24x
differentiate wrt x
0.956 s0.031 s
30x
compute f (6, y, z)
3.760 s0.175 s21x
reorder exponents [d, x, y, z] → [d, y, z, x], sort, Horner form
multiply by a term
1.190 s0.066 s18x
search for subexpression
0.040 s0.002 s20x
linear search to degree 100, O(n) for names, O(log t) for terms
set of indeterminates
0.060 s0.000 s→ O(1)
leading coefficient in x
0.058 s0.005 s11x
monomial division avoids comparisons, stop early using degree
extract terms of f
0.634 s2.420 s0.26x
create old structure, cached until the next garbage collection
replace variable
1.160 s0.076 s15x
add exponents of like variables, sort and combine like terms
Taylor series to O(x50 )
0.668 s0.055 s12x
reorder exponents [d, x, y, z] → [x, d, y, z], sort to collect in x
type check0.029 s0.000 s→ O(n)
type check the variables, other types are optimized as well

For some algorithms we need to sort polynomials after modifying themonomial data. We sort by treating the monomials as unsigned 64-bit integers and using an inplace MSD American flag radix sort [15]. For example,
the coeffs(f,x), eval(f,x=6) and taylor(f,x,n) commands all need thecoefficients of f in x. Suppose f is a polynomial in {w, x, y}. For each monomial wi xj y k in f , encoded as dijk where d = i + j + k is its degree, with

POLY : A new polynomial data structure for Maple 179

a constant number of masks and bit operations (seven suffices) we move j,
the degree of x to the front to obtain the jdik . We sort this modified data
to group terms by power of x with ties broken by the monomial ordering on{w, y}. From that we extract the coefficients in x, already sorted, in lineartime.
The biggest improvement we have seen for a Maple library command is thecollect command, which is used to write a polynomial in recursive form. For
example, if f = xy 3 + x2 y − x2 z + xyz − 2 then collect(f,x) will rewrite fas (y − z)x2 + (y 3 + yz)x − 2. The Maple code for the collect command uses
the series(f,x,3) command to implement this. Since the series command
is 8x faster, we did not expect collect to be 6.196/0.111 = 56 times faster.
Below is a profile showing that most of the time in Maple 16 was not inseries at all, but rather in the frontend and indets commands.
Table 3 Profiles for executing collect(f,x) in Maple 16 and Maple 17.
function
frontend
indetsseriescollect
total:

calls
1
2
1
1
5

Maple 16
time time%
3.672
59.26
1.617
26.10
0.747
12.06
0.160
2.58
6.196 100.00

Maple 17time time%
0.004
3.60
0.000
0.00
0.094
84.68
0.013
11.71
0.111 100.00

In Maple 17, the cost of frontend and indets are now negligible since
they no longer need to descend into the sum-of-products dag. In Maple 16,
frontend recursively descends the sum-of-products dag looking for objects
like x1/2 , sin(x), 2n , etc. But our polynomial has none; it only has variables
x, y and z in it.

3.1 Unpacking
The one case where we lose is when we must unpack the POLY dag andconvert to the old data structure. The Maple command op(f) constructs asequence of all the terms of f . Other Maple commands which effectively do
the same thing include the common Maple programming idioms:
map(g,f)
for t in f doindets(f,t)

apply the function g to each term of f
iterate over the terms of f
extract all subexpressions in f of type tEach term, e.g. 8xy 2 , is stored as SUM ↑ P 8 where P is the monomial stored
as PROD x 1 y 2 . Thus Maple 17 must build a SUM and a PROD for eachterm in f whereas Maple 16 only builds a SUM and the PROD already exists.

10Michael Monagan and Roman Pearce

Theoretically, Maple 17 is O(nt) compared with O(t) for Maple 16. We have
tried to improve the speed of unpacking by creating the PROD objects in
simplified form, but the slowdown for op(f) remains a factor of 4.
However, we observed a tendency of Maple library code to either frequently
unpack POLY or rarely unpack, often in the course of checking high leveltypes. Alongside other internal caches in the Maple kernel (e.g. for subs andindets) we added a cache for unpacked POLY dags. This cache is clearedout on every garbage collection, so its practical effect is restricted to small
polynomials that would be repeatedly unpacked.
For indets(f,t) we can avoid unpacking in many cases by detecting types
that do not appear in POLY or appear only in the variables. Table 4 shows
the most common types in indets(f,t) when f is a POLY when the entire
Maple library test suite is run. For products or powers we must unpack theterms of f to create the result, but in the top ten cases (and many others)
that is avoided.
Table 4 Calls to indets(f, t) for f a POLY in the Maple library test suite.
type tnamenonrealfloat
function
Or(RootOf, radical){name, function}

number11937973
708148669307776678146
1863699
1861368

type t{rtable, table}
specfunc(anything, RootOf)
radicalindexed‘^‘Or(‘+‘,‘*‘,‘^‘)

number1509366
1429737
11015391089504
1047368828257

3.2 Repacking
A number of the Maple kernel operations require us to repack monomials. For
example, when adding x2 + y 2 + z 2 and x2 + y 2 − z 2 , the result 2x2 + 2y 2 does
not have the variable z. When Maple simplifies the result it must detect thatz is missing and repack the polynomial into Z[x, y]. Repacking also occurs in
the coeff, coeffs, eval, lcoeff, and taylor commands that remove one
or more variables, or expand and divide, which convert polynomials to acommon ring. Critically, these operations do not permute the variables; they
insert or remove exponent fields and change the sizes of the remaining fields.
In Maple 17 this is coded as a straightforward O(n) loop, but we digresson the topic because new Intel microprocessors (codenamed Haswell, see [13])have added two new instructions with exciting possibilities. The PEXT instruction is short for parallel extract, and it flushes masked bits to the bottomof the word.

POLY : A new polynomial data structure for Maple 17

11

monomial x3 y 5 z 7 00001111 00000011 00000101 00000111
mask for {x, z} 00000000 11111111 00000000 11111111
result from PEXT 00000000 00000000 00000011 00000111
Its inverse is PDEP, short for parallel deposit, which distributes bits to various
locations in a word, starting from the lowest bit. Both operations are O(1).
input data 00000000 00000000 00000011 00000111
mask of locations 00000000 11111111 00000000 11111111
result from PDEP 00000000 00000011 00000000 00000111
The parallel design for these operations was originally proposed by Hilewitz
and Lee in [12]. With these instructions, we can repack monomials without
any branches or loops, e.g. to convert from Z[x, z] to Z[x, y, z] we would do:
3 × 21 bits
d i kextract−→

4 × 16 bits
0 d i kdeposit−→4 × 16 bits
d i 0 kThis would be needed, for example, to multiply f (x, z) and g(x, y) inZ[x, y, z]. In fact Maple 17 does use the PEXT operation to reorder multiple variables for coeffs, but we coded this in C using the algorithm inHacker’s Delight [21]. After precomputation it uses 24 bitwise instructions
per word, and it is called twice per term by coeffs to reorder the exponents.
We attempted to use a C routine for PDEP to repack monomials but it was
hard to get a gain. Nevertheless this should be our approach as soon as there
is widespread hardware support.

3.3 Hashing and Simplification
When an algebraic expression is created in Maple, it is simplified recursively
by the kernel. For example, consider the polynomial f = 9xy 3 z − 4y 3 z 2 −6xy 2 z − 8x3 − 5 from Figure 1. Each object in f that is not a small integer
is simplified and hashed to see if it already exists in memory. In Maple 16,
this is first done for the variables x, y, z then for the monomials (PRODs in
Figure 1) xy 3 z, y 3 z 2 , xy 2 z and x3 , and finally for the whole expression f .
These objects are hashed and stored in the internal simpl table, which maps
each object to a unique copy in memory. This feature of Maple allows it toidentify equal objects by address.
In Maple 17 (see Figure 3), because monomials are encoded immediately
as machine integers, they are not stored in the simpl table. Rather, only
the POLY object, the SEQ object (the sequence of variables), and multiprecisioncoefficients are stored in the simpl table.
What is gained from not having to create, simplify, and hash each monomial as a PROD object? The following benchmark gives us a clue. Consider> f := expand((1 + s + t + u + v + w + x + y + z)16 ):12Michael Monagan and Roman Pearce

which creates a polynomial f in 8 variables with 735471 terms. By first issuing
the command sdmp:-info(1): we can obtain profiling information from the
C library that computes this result. Table 5 shows that most of the time isspent simplifying the result in Maple 16, whereas in Maple 17 this time isfairly small. The C routine itself is identical with only a tiny difference dueto the compiler.
Table 5 Real time in seconds for computing and simplifying a large power.
Maple 16
Maple 17

expand power allocate dag simplify dag total time
0.133 s0.080 s1.180 s1.420 s0.128 s0.000 s0.010 s0.139 sThe C library uses the same monomial representation as Maple 17, so it
can copy the term data to a new POLY object. Maple 16 has to allocate thesum of products structure, which is almost as expensive as the computation.
Maple 17 simplifies the POLY by checking that its terms and variables are
sorted and all variables have a non-zero exponent. This takes 0.01 seconds or
7.2% of the time. Maple 16 must do considerably more work to simplify thesum of products dag. For each PROD, it checks that the variables are distinct(they are) using an O(n2 ) loop, then it hashes the PROD and inserts it into
the simpl table. Then it has to sort the SUM to check that all the PRODs
are distinct (they are), because hashing has destroyed any previous order.
Finally, it hashes the SUM.
A hidden cost is that the code to simplify SUMs and PRODs is quite expensive. It handles special objects like infinity and undefined, complex andfloating point arithmetic, operator overloading, binary relations like < or =because equations can be added or scaled, matrix arithmetic, etc. These routines implement much of the expressive power of the Maple language, and thisis not free. To simplify POLY, we have the luxury of analyzing the variables
and calling algorithms that work in restricted domains.

4 Benchmarks
What is the impact of the POLY data structure on Maple’s overall performance? This was difficult to predict in advance. One goal was to reduce
sequential overhead in polynomial algorithms so that parallel speedup in multiplication and division (see [18, 19]) would speed up the Maple library. Was
that achieved? To this end we developed two benchmarks; expanding determinants of polynomial matrices and factoring multivariate polynomials. Bothare higher level algorithms.

POLY : A new polynomial data structure for Maple 17

13

4.1 Determinant Benchmark
Our first benchmark computes the determinant of the n × n symmetric
Toeplitz matrix A for 6 ≤ n ≤ 11. This is a matrix with n variables{x1 , . . . , xn } with xi along the ith diagonal and ith subdiagonal. To compute det(A) we use our own implementations of the Bareiss algorithm (see[2]), which we provide in Appendix A. At the k th elimination step, ignoring
pivoting, the Bareiss algorithm computes
Ai,j :=

Ak,k Ai,j − Ai,k Ak,jfor i = k + 1, . . . , n and j = k + 1, . . . , n (1)Ak−1,k−1

where the division is exact. At the end of the algorithm An,n = ± det(A).
Thus the Bariess algorithm does a sequence of O(n3 ) polynomial multiplications and divisions that grow in size, with the largest one occurring in thelast step when k = n − 1.
In Table 6 below, #det is the number of terms in the determinant which has
degree n, and #num is the number of terms in An−1,n−1 An,n −An,n−1 An−1,nwhich has degree 2n − 2 and is much larger than det(A).

Table 6 Real time in seconds to compute det(A) using the Bareiss algorithm.
n #det #num
6
120
575
7
427
3277
8 1628 21016
9 6090 128530
10 23797 813638
11 90296 5060172

Maple 13
Maple 16
Maple 17
Magma 2.17
1 core 1 core 4 cores 1 core 4 cores
1 core0.015 0.008 0.009 0.002 0.002
0.000
0.105 0.030 0.030 0.010 0.010
0.020
1.123 0.181 0.169 0.047 0.037
0.200
19.176 1.450 1.290 0.482 0.294
2.870
445.611 14.830 12.240 5.680 2.940
77.020−151.200 94.340 78.500 26.890
2098.790

We used a quad core Intel Core i7 920 2.66 GHz CPU running 64-bit
Linux. Timings are real times in seconds. With four cores we achieve a factorof 3 to 4 speedup over Maple 16, which is large. That gain is entirely fromreducing the overhead of Maple data structures; Maple 16 and Maple 17
use the same C library routines for polynomial arithmetic. The reduction ofoverhead increases parallel speedup from 1.6x to 2.59x over Maple 16. For
comparison we include times for Maple 13 (which does not use our C library)and Magma 2.17.
By default, Maple and Magma do not use the Bareiss algorithm to compute these determinants. Instead, they use the method of minor expansion of
Gentleman and Johnson [7]. Recall that given an n × n matrix A
det(A) =n
X
i=1(−1)n+1 Ai,1 det(M (1, i))

(2)

14

Michael Monagan and Roman Pearce

where M (1, i) is the n − 1 by n − 1 matrix obtained from A by deletingcolumn 1 and row i. Applied naively, this formula recomputes the determinants of sub-matrices many times. Gentleman and Johnson avoided that bycomputing from the bottom up; they compute all 2 × 2 determinants thenall3 × 3 determinants and so on. This is still exponential
 computes
Pnin n.nItnndeterminantsofk×k
sub-matrices
for
a
total
of
k=1 k = 2 − 1kdeterminants.
For our Toeplitz matrices, the multiplications in (2) are of the form variable
times polynomial. For example, to multiply f = 9xy 3 z−4y 3 z 2 −6xy 2 z−8x3 −5
by y, we add the monomial representation for y to each monomial in f , namely
y = 1010 = 248 + 216 to the encodings of 5131 , 5032 , 4121 , 3300 , 0000 .
Notice how the result remains sorted in the monomial ordering. The additionsin (2) are n-ary, for which we copy terms to a new POLY and apply radix
sort. The improvement shown in Table 7 was huge. It surprised us.
Table 7 Real times in seconds for determinants using minor expansion.
n
6
7
8
9
10
11
12
13

#det
120
427
1628
6090
23797
90296
3507261338076

Maple 16
0.002
0.010
0.049
0.305
1.991
19.370274.9902024.370

Maple 17 Magma 2.17
0.002
0.001
0.004
0.003
0.013
0.019
0.047
0.116
0.252
0.770
1.322
6.210
6.737
44.50037.570337.770

4.2 Factorization benchmark
Our second benchmark is multivariate factorization. For perspective we include timings for Magma [3], Mathematica, Maxima [14], Sage [20], Singular[10], and Trip [5] which is a computer algebra system for celestial mechanics.
Table 8 reports the real times for multiplication, division, and factorization
on a hyperthreaded quad core Intel Core i7 920 2.66 GHz running 64-bit
Linux. For each timing we report the median of three test runs. Maple 16
and 17 start up to four threads depending on the size of each multiplication
or division. The factorization routine is sequential Maple code, which gains
parallelism from our multiplication and division routines.
For Mathematica 9 we timed the internal functions Algebra‘IPExpand
and Algebra‘IPExactQuotient for multiplication and division. The additional time for the top level Expand was small, and we found no suitablecommand for exact division at the top level. The Factor command did not
make use of parallelism.

Maple 13
multiplyp1 := f1 · (f1 + 1)
p2 := f2 · (f2 + 1)
p3 := f3 · (f3 + 1)
p4 := f4 · (f4 + 1)
p5 := f5 · g5
p6 := f6 · g6
divideq1 := p1 /f1
q2 := p2 /f2
q3 := p3 /f3
q4 := p4 /f4
q5 := p5 /f5
q6 := p6 /f6factorp1 : 12341 terms
p2 : 12341 terms
p3 : 38711 terms
p4 : 135751 terms
p5 : 12552 terms
p6 : 417311 terms

Maple 16
1 core 4 cores

Maple 17Magma
1 core 4 cores2.19-11.561
1.544
26.50198.351
13.666
11.486

0.063
0.063
0.521
2.180
1.588
0.772

0.030
0.032
0.171
0.649
0.384
0.628

0.041
0.041
0.403
1.814
0.153
0.204

0.012
0.012
0.102
0.416
0.154
0.082

0.33 1.09 0.57 0.585
0.33 1.08 0.58 0.608
3.99 10.86 6.77 6.595
13.70 61.77 30.99 31.806
13.24 32.73 18.22 17.776
0.89 3.04 2.75 1.787

0.56
5.92
0.57
5.86
24.25 59.32
5.83 306.10
15.83 196.28
2.52 22.45

0.120
0.115
0.855
5.732
1.526

0.011
0.031
0.094
0.488
0.177
0.064

0.008
0.015
0.080
0.388
0.145
0.047

1.451
1.435
25.05492.8675.570
10.421

0.065
0.065
0.524
2.253
1.636
0.769

0.033
0.033
0.184
0.736
0.417
0.627

0.042
0.042
0.411
1.842
1.445
0.215

0.015
0.015
0.117
0.483
0.333
0.095

0.36 1.30 0.40 0.183
0.36 1.30 0.39 0.183
4.14 17.81 3.64 4.737
18.54 80.39 14.96 11.420
12.48 50.16 10.80 10.478
7.90 4.87 1.84 1.484

0.55
5.00
0.56
5.04
24.63 47.01
5.93 228.83
16.14 35.94
2.69 30.28

0.197
0.194
1.510
5.662
2.789

0.238
0.225
2.340
9.540
8.827
1.181

0.191
0.208
1.776
7.374
5.670
1.016

31.330275.508360.8622856.388
302.453
1359.473

2.792
3.240
16.71459.00926.43551.7022.658 0.790
3.071 0.991
14.110 6.927
46.151 24.345
16.152 12.131
48.808 8.295

0.650
0.850
4.399
12.733
6.800
6.330

6.51
7.09
119.32320.04105.55369.1218.478112.243276.161951.725935.149

NA
NA
NA
NA
NA
NA

NA
NA
NA
NA
NA
NA

f1 = (1 + x + y + z)20 + 1
1771 terms

f2 = (1 + x2 + y 2 + z 2 )20 + 1
1771 terms

f5 = (1 + x)20 (1 + y)20 (1 + z)20 + 1
g5 = (1 − x)20 (1 − y)20 (1 − z)20 + 1
9261 terms

Sage
SingularMaxima5.8 3-1-4 3-1-6
5.25.0

1.51
1.58
18.14
68.32
14.63
37.56

2.01 0.853
2.10 0.933
12.48 10.074
61.85 39.353
13.83 9.604
42.08 20.603

4.54
4.97
163.0644.94
1046.67
155.49f3 = (1 + x + y + z)30 + 1
5456 terms

Mathematica
8.0
9.0

11.70
67.83
166.10604.74563.54290.07

Trip 1.2.26(RS) (RD)POLY : A new polynomial data structure for Maple 17

Table 8 Real times in seconds for polynomial multiplication, division and factorization.

f4 = (1 + x + y + z + t)20 + 1
10626 terms

f6 = (1 + u2 + v + w2 + x − y)10 + 1
g6 = (1 + u + v 2 + w + x2 − y)10 + 1
3003 terms
15

16

Michael Monagan and Roman Pearce

We report two timings for Trip: the (RS) time is for the optimized recursive sparse polynomial data structure POLYV, while the (RD) time is the
optimized recursive dense data structure POLPV. Both use Trip’s parallel
routines (see [6]) with 8 threads and rational arithmetic, including a fast
representation for small machine integers similar to Maple’s.
There are some anomalies in Table 8. Maple’s times for division are close
to those for multiplication, except on Problem 5 where Maple 17 uses a
sequential dense method to multiply. Singular’s timings for division are often
twice as fast as multiplication. This is because Singular multiplies in the
distributed representation and divides in a recursive representation. Trip’s
times for division are much slower than those for multiplication. This is partly
because division in Trip is not parallelized.
In comparing the timings for factoring p1 and p2 we see that factoring p2
is much slower in Maple 13 and Mathematica, but not in Maple 16 or 17,
Magma, or Singular. The fast systems apply the substitution p2 (x2 = u, y 2 =v, z 2 = w) to reduce the degree of the input polynomial before factoring it.
This halves the number of Hensel lifting steps in each variable.
We note that Singular’s timings for factorization have improved enormously since version 3-1-0. Times for version 3-1-0 on the first four factorizations were 12.28, 23.67, 97.10, 404.86 seconds. The factorization code in 3-1-4
was changed by Michael Lee to use a recursive representation for polynomials
instead of the default distributed representation shown in Figure 2.
Our first improvement from Maple 13 to Maple 16 was due to our C libraryfor polynomial multiplication and division described in [16, 17, 18, 19]. The
speedup in multiplication and division produces a speedup in multivariate
factorization, because most of the time is spent in “Hensel lifting” which
consists of many multiplications and some exact divisions.
Our second improvement was to parallelize the algorithms for multiplication and division. In many cases we obtain superlinear speedup in our C
library but the top level expand and divide have lower speedup, because of
the extra time to import and export Maple data structures. For higher level
algorithms such as factor, parallel speedup is further reduced by the needto perform many small operations in sequence. The cost of degree, indets,
and type(f,polynom) also reduce parallelism in higher level code.
With the introduction of the POLY dag in Maple 17, we have substantially reduced data structure overhead and the cost of almost all supporting
routines. Table 9 shows the improvements to parallel speedup which come on
top of the large gains achieved for sequential time. The speedup for expand
and divide is now much closer to our C library, and the speedup for factor,
while modest, is respectable for a sequential algorithm.
Table 10 shows the speedup of Maple 17 over Maple 16, on 4 cores versus
1. Notice how the gain is larger in parallel. This is just a pleasant consequenceof Amdahl’s Law when you reduce sequential overhead in parallel algorithms.
The sole exception (Problem 5 expand) uses a sequential dense method in
Maple 17.

POLY : A new polynomial data structure for Maple 1717

Table 9 Parallel speedup (1 core)/(4 cores) in Maple 17 versus Maple 16.
expanddividefactor

Maple 17
3.41 3.41 3.95 4.36 1.00 2.49
2.80 2.80 3.51 3.81 4.34 2.26
1.21 1.16 1.57 1.91 1.78 1.31

Maple 16
2.10 1.97 3.04 3.36 4.13 1.23
1.97 1.97 2.85 3.06 3.92 1.23
1.05 1.05 1.18 1.28 1.63 1.06

Table 10 Observed speedup (Maple 17)/(Maple 16) on 4 cores versus 1 core.
expanddividefactor

4 cores
2.50 2.66 1.67 1.56 2.49 7.65
2.20 2.20 1.57 1.52 1.25 6.60
4.09 3.61 3.20 3.62 2.37 7.71

1 core1.53 1.53 1.29 1.20 10.38 3.78
1.54 1.54 1.27 1.22 1.13 3.57
3.53 3.27 2.41 2.42 2.18 6.32

With the POLY dag in Maple 17, the time for factorization on one core has
been reduced by more than 50%, but the parallel speedups are even greater:
The savings are entirely sequential time, as can be seen by subtracting theparallel times from the sequential times. The reduction in overhead improves
parallel speedup significantly, from 1.28x to 1.91x in the case of factoring p4 .
To see where the improvements in the factorization have come from wehave profiled the main parts of the factorization code. The profile (see Table 11) shows the %age of the time in the main parts of the factorization
algorithm for Maple 16 and Maple 17. The data under improved coeftayl
includes a further algorithmic improvement. The data shows we have eliminated 0.599−0.377 = 0.222s of overhead from the polynomial multiplications(see row expand) or 37%. The biggest speedup is division (see row divide).
This is because the divisions are mostly trial divisions which fail quickly. Insuch cases almost all the time is in conversion which is wasted.
Table 11 profile for factor(p1); (1 core).
Maple 16
Maple 17
improved coeftayl
functiontime time% time time% time
time%
coeftayl
1.086s41.06 0.310s
28.21 0.095s
12.03
expand
0.506s19.13 0.263s
23.93 0.255s
32.28
diophant
0.424s
16.03 0.403s
34.94 0.299s
37.85
divide0.256s
9.68 0.034s
3.09 0.035s
4.43
factor
0.201s7.60 0.011s
1.00 0.010s
1.27
factor/hensel
0.127s4.80 0.064s
5.82 0.063s
7.97
factor/unifactor 0.045s
1.70 0.033s
3.00 0.033s
4.18
total:
2.645s 100.00% 1.099s 100.00% 0.790s
100.00%The biggest absolute gain is for the routine coeftayl(f,x-a,k) which computes the coefficient of f in (x − a)k . This computation is not done by expanding f as a Taylor series about x = a, but instead, by using the formulag(x = a)/k! where g = ddf
k x , the k’th derivative of f . Referring back to Table

18

Michael Monagan and Roman Pearce

2, we can see that the speedup is due to the improvement of differentiationand polynomial evaluation. We also tried the following formula to computePdegx fthe coefficient:coeff(f, xi )ai ki . We can see that this is 3× fasteri=kagain (see improved coeftayl). The total real time is reduced from 2.59s to
1.07s to 0.790s.

5 Conclusion
Maple, Mathematica, Magma and Singular all use a distributed representation for multivariate polynomials. Maple’s sum-of-products data structureand Singular’s linked list data structure are illustrated in Figures 1 and 2
in the introduction. We ask the the reader take another good look at them.
Mathematica’s data structure is similar to Maple’s and Magma’s data structure is similar to Singular’s. These data structures, which were designed in the
1980s when memory access was constant time, will not yield high-performanceon todays computers because memory access is not sequential.
One way to speed up polynomial multiplication, division, or factorization
would be to convert the input to a more suitable data structure, compute theresult, then convert back. This is what we did in [16] for Maple 14 for polynomial multiplication and division. Singular 3-1-4 does this for polynomial
division and factorization. It switches to using a recursive representation fordivision and factorization. However, conversion overhead will limit parallel
speedup. Amdahl’s law states that if the sequential proportion of a task is S
then parallel speedup on N cores is limited to
speedup ≤1.
S + (1 − S)/N

When S is large (50% or more say), then we cannot get good parallel speedup.
What we have done in this work for Maple is to make our POLY data
structure the default data structure in Maple. The POLY data structureis used when all monomials in a polynomial can be packed into a single
word. This enabled us to eliminate conversion overhead in multiplication anddivision. The data in Table 6 shows improved parallel speedup for polynomial
multiplication and division. We also implemented highly efficient algorithms
for many Maple kernel operations for POLY. The data in Table 7 shows aspeedup of a factor of 50 over Maple 16 for a routine polynomial determinantcomputation. The data in Table 8 shows speedups of factors of between 2
and 3 for large multivariate polynomial factorizations which is a huge gain.
Although not reported here, we also find speedups of a factor of 2 for large
multivariate polynomial gcd computations.
The cost incurred is mainly in code complexity. We must manage twodata structures for polynomials, one where the coefficients are integers andPOLY : A new polynomial data structure for Maple 17

19

the monomials can be packed into a single machine word, and one, Maple’s
sum-of-products data structure, which does not have these restrictions. A
substantial programming effort was required to support the new data structure in the Maple kernel. The gains suggest this is worthwhile.
In closing, the reader may have wondered why we only use one word ofmemory to encode monomials, and not two, or more? For if we use twowords, we could encode polynomials in twice as many variables or of muchhigher degree. With 128 bits, one will cover almost all applications. We would
like to see how far 64-bits takes us before considering such an extension.
For supporting two word exponents potentially doubles the amount of code.
Another desirable extension is to allow the coefficients in the POLY dag tobe fractions or floating point numbers as well as integers.

References
1. O. Bachmann and H. Schönemann. Monomial representations for Grobner bases computations. Proceedings of ISSAC ’98, pp. 309–316, 1998.
2. E. Bariess, 1968. Sylvester’s Identity and Multistep Integer-Preserving Gaussian Elimination. Mathematics of computation 22 (102): 565−578.
3. Bosma, W., Cannon, J., Playoust, C., 1997. The Magma Algebra System I: The User Language. J. Symb. Cmpt. 24(3-4), 235–265. See also
http://magma.maths.usyd.edu.au/magma
4. Buchberger, B. Private Communication, May 2013.
5. Gastineau, M., Laskar, J., 2006. Development of TRIP: Fast Sparse Multivariate Polynomial Multiplication Using Burst Tries. Proceedings of ICCS 2006, Springer LNCS
3992, pp. 446–453.
6. Gastineau, M., 2010. Parallel operations of sparse polynomials on multicores - I. Multiplication and Poisson bracket. Proceedings of PASCO ’2010, ACM Press, pp. 44–52,
2010.
7. Gentleman, W.M., Johnson, S.C. Analysis of Algorithms, A Case Study: Determinantsof Matrices with Polynomial Entries. ACM Trans. on Math. Soft., 2(3), pp. 232–241,
September 1976.
8. Granlund, T., 2008. The GNU Multiple Precision Arithmetic Library, version 4.2.2.
http://www.gmplib.org/
9. Grayson, Daniel R. Stillman, Michael E., Macaulay2, a software system for researchin algebraic geometry. Available at http://www.math.uiuc.edu/Macaulay2/
10. Greuel, G.-M., Pfister, G., Schönemann, H., 2005. Singular 3.0: A Computer Algebra
System for Polynomial Computations. Centre for Computer Algebra, University of
Kaiserslautern. http://www.singular.uni-kl.de
11. Hall, A.D. Jr., The ALTRAN System for Rational Function Manipulation – A Survey.
Communications of the ACM, 14, 517–521, ACM Press, 1971.
12. Hilewitz, Y., Lee, R.B., 2006. Fast Bit Compression and Expansion with Parallel Extract and Parallel Deposit Instructions. Proceedings of ASAP ’06, IEEE, pp. 65–72,
2006.
13. Intel Corporation. Advanced Vector Extensions Programming Reference. June 2011.
http://software.intel.com/sites/default/files/m/8/a/1/8/4/36945-319433-011.pdf
14. The Maxima computer algebra system. http://maxima.sourceforge.net/
15. Peter M. McIlroy, Keith Bostic, and M. Douglas McIlroy. Engineering Radix Sort,
Computing Systems, 6(1): 5–27, 1993.

20

Michael Monagan and Roman Pearce

16. Monagan, M., Pearce, R., 2010. Sparse Polynomial Multiplication and Division in
Maple 14. Communications in Computer Algebra, 44:4, 205–209, December 2010.
17. Monagan, M., Pearce, R., 2011. Sparse Polynomial Division using Heaps. J. Symb.
Cmpt. 46(7):807–822, 2011.
18. M. Monagan, R. Pearce., 2009. Parallel Sparse Polynomial Multiplication Using Heaps.
Proceedings of of ISSAC 2009, ACM Press, pp. 295–315.
19. M. Monagan, R. Pearce., 2010. Parallel Sparse Polynomial Division Using Heaps. Proc.
of PASCO 2010, ACM Press, pp. 105–111.
20. William A. Stein et. al. Sage Mathematics Software (Version 5.0), The Sage Development Team, 2012, http://www.sagemath.org.
21. Warren, Henry S. Hacker’s Delight. Addison-Wesley, 2003.

Appendix A
Maple code (no pivoting) for the Bareiss algorithm.
ffge := proc(A,n) local d,i,j,k,t;
d := 1;for k to n-1 dofor i from k+1 to n do
for j from k+1 to n dot := expand(A[k,k]*A[i,j]-A[i,k]*A[k,j]);divide(t, d, evaln(A[i,j]));od;
A[i,k] := 0;od;d := A[k,k];od:
A[n,n];
end;
n := 8;
T := linalg[toeplitz]([seq(x[i],i=1..n)]);A := array(1..n,1..n):for i to n do for j to n do A[i,j] := T[i,j] od od:
det := CodeTools[Usage]( ffge(A,n) ):Magma code for the Bareiss algorithm.
Z := IntegerRing();
P<x,y,z,u,v,w,p,q,r,s,t,a> := PolynomialRing(Z,12);
X := [x,y,z,u,v,w,p,q,r,s,t,a];n := 8;
A := Matrix(P,n,n,[0 : i in [1..n^2]]);for i in [1..n] dofor j in [1..n] do
A[i,j] := X[AbsoluteValue(j-i)+1];

POLY : A new polynomial data structure for Maple 17

21

end for;end for;d := 1;time for k in [1..n-1] dofor i in [k+1..n] dofor j in [k+1..n] dot := A[k,k]*A[i,j]-A[i,k]*A[k,j];A[i,j] := ExactQuotient(t,d);end for;end for;d := A[k,k];
end for;det := A[n,n];Appendix B
Maple code for timing benchmarks.
f := expand( (1+x+y+z)^20 )+1:
p := CodeTools[Usage]( expand( f*(f+1) ) ):CodeTools[Usage]( divide(p,f,’q’) );CodeTools[Usage]( factor(p) ):Magma code for timing benchmarks.
Z := IntegerRing();
P<x,y,z> := PolynomialRing(Z,3);
f := (1+x+y+z)^20+1;
g := f+1;
time h := f*g;time q := ExactQuotient(h,f);time ff := Factorization(h);Mathematica code for timing benchmarks. Note, for the fifth benchmark, weused Cancel[p/f] which was much faster than PolynomialQuotient.
f = Expand[(1+x+y+z)^20]+1;
AbsoluteTiming[p = Expand[f*(f+1)];]
AbsoluteTiming[q = PolynomialQuotient[p,f,x];]
AbsoluteTiming[h = Factor[p];]Maxima code for timing benchmarks.

22

Michael Monagan and Roman Pearce

showtime : true;f : rat( (1+x+y+z)^20 ) +1 $h : f*(f+1)$
qr : divide( h, f )$f : factor( h )$Sage code for timing benchmarks
Q = RationalField()
P.<x,y,z> = PolynomialRing(Q,3,order=’deglex’)f = (1+x+y+z)^20+1%time p = f*(f+1)%time q,r = p.quo_rem(f)%time h = factor(p)Singular code for timing benchmarks.
ring R=0,(x,y,z),lp;
poly f = (1+x+y+z)^20+1;
poly g = f+1;
int TIMER;
TIMER = timer; poly p = f*g; timer-TIMER;
TIMER = timer; poly q = p/f; timer-TIMER;
TIMER = timer; list L = factorize(p); timer-TIMER;
Trip code for timing benchmarks. POLYV means recursive sparse, POLPV
means recursive dense.
reset; _cpu=8$ _mode=POLYV$ _modenum=NUMDBL$
f=(1+x+y+z)^20+1$ g=f+1$ p = 0$
time_s; p = f*g$ time_t(usertime, realtime); realtime;
time_s$ div(p,f,q,r)$ time_t(ctime,rtime)$ rtime;