POLY: A New Polynomial Data Structure for Maple 17

Michael Monagan

POLY: A New Polynomial Data Structure for Maple 17

2014, Computer Mathematics

Abstract

We demonstrate how a new data structure for sparse distributed polynomials in the Maple kernel significantly accelerates several key Maple library routines. The POLY data structure and its associated kernel operations (degree, coeff, subs, has, diff, eval, ...) are programmed for high scalability with very low overhead. This enables polynomial to have tens of millions of terms, increases parallel speedup in existing routines and dramatically improves the performance of high level Maple library routines.

Figures (9)

Figure 1: Maple’s sum-of-products representation has irregular Maple dags for each term. Figure 1 below shows the default polynomial data structure in Maple 16 and all previous versions for the polynomial f = 9xy?z— 4y?z? — 6axy?z— 8x3 —5. It is a “sum-of-products” where each term has a separate Maple object, a PROD, to represent the monomial. To compute the degree of f, a coefficient in x, test for a subexpression, or do almost anything else, the Maple kernel must recursively descend through multiple levels of dags. This involves extensive branching and random memory access, both of which are slow, and will prevent Maple from achieving high-performance on modern CPUs.

Figure 2: The new packed distributed representation. Figure 2 below shows our new data structure for sparse distributed polynomials. The first word is a header word. The second word is a pointer 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 together with the total degree in a single machine word. E.g. for xy?z? we store the values (6,1,2,3) as 6-248 4 992 4.9.91 4 3 on a 64-bit machine. The terms are sorted into graded lex order by comparing the monomials as unsigned integers. This gives a canonical representation for the polynomial. A small integer coefficient —2°? < x < 2° is encoded as 2% + 1 so that the right most bit is 1. Large integer coefficients are encoded as a pointer (with right most bit 0) to a GMP integer ((7]).

Table 1: Improvements for Maple kernel operations. The new representation has allowed us to write many 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. n the new data structure, we can often avoid doing expensive operations on all of the terms, or we can do them much more efficiently. In Table 1 we measure the speedup on a polynomial with one million terms in three variables, constructed as f := expand(mul(randpoly(i, degree = 100, dense),i = [x,y, z])): The cost or evaluation is added to the other commands if you are using Maple interactively. Note, none of these operations in either Maple is parallelized.

Table 2: Profile for executing collect (f,x) in Maple 16 and our new Maple. The biggest improvement we have seen for a Maple library command is the collect command. It is used to write a polynomial in recursive form. For example, if f = xy? +a?y—2?z+ayz—2, the command collect (f,x) writes f as (y—z)a?+(y?+yz)a—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 in our new Maple, we were not expecting collect to be 6.616/0.123 = 54 times faster. Here is a profile showing that most of the time in Maple 16 was not in the series command, but in the frontend and indets commands. In our new Maple, the cost of frontend and indets are now negligible since they no longer need to descend into the sum-of-products dag. In the new Maple, they only need to look at the variables which costs O(n). Why were frontend and indets so expensive? They need to search the sum-of-products dag to see if there are any indeterminates which are not variables. They are looking for objects like x!/?, sin(2), 2”, etc. But our polynomial has none; it only has variables x,y and z in it. In order to do this they pick apart each product and each power, e.g., given z°yz* they recursively constructs x? and z4 as new objects before descending them to see 2, 3, z, 4.

Table 3: timings (real times in seconds) for determinants using the Bareiss algorithm.

Table 4: timings (real times in seconds) for determinants using minor expansion all 2 x 2 then all 3 x 3, etc. This is still exponential in n; it computes (,) determinants of sub-matrices ot size k by k for a total of >y_, (;) = 2" — 1 determinants. For our Toeplitz matrices, the For example, to multiply nomial. representation for y, namely 1010 polynomial multiplications in (2) are all of the form variable x poly- f = 9xyz — 4y?z? — 6ry?z — 82° — 5 by y, we add the monomial = 248 + 2! to each monomial in f, namely the integer encodings of 5131 |, 5032 > 4121 |, 3300 2 0000 |. Notice that result will already be sorted in the graded ordering. The polynomia additions and subtractions in (2) are also done in linear time by a simple merge. The improvement shown in Table 4 below is huge. It surprised us.

There are some anomalies in Table 5. Maple’s timings for division on 1 core and 4 cores are very close to those for multiplication. However Singular’s division timings for p3/f3 and p4/f4 are more than twice as fast as the time for multiplication. This is because Singular uses a distributed data structure to represent polynomials. It multiplies in this distributed representation but for division, switches to using a recursive representation. On the other hand, Trip’s timings for division are much slower than for multiplication. This is partly because division in Trip 1.2 has not been parallelized. In our second benchmark we see a large gain in performance on polynomial factorization. To provide some perspective, we include timings for Magma [3], Singular [8], Mathematica, and Trip [4], a computer algebra system for celestial mechanics.

a a a i SS YER eM’ To see where the SpTOTETGENIS in the factorization have come from we have profiled the main parts of the factorization code. The profile (see Table 6) shows the %age of the time in the main parts of the factorization algorithm for Maple 16 and our new Maple. 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. In such cases almost all the time is in conversion which is wasted.

Log In

POLY: A New Polynomial Data Structure for Maple 17

Sign up for access to the world's latest research

Abstract

Related papers

Related topics

Related papers