Our first article introduced Maple's polynomial data structure and explained how Maple spends a lot of time working with monomials. To multiply polynomials having n and m terms, Maple must construct, simplify, hash, and sort all nm pairwise products to determine what monomials are equal. This work is performed even if the result has far fewer than nm terms, making it a rather inefficient way to multiply large multivariate polynomials. This article describes a new data structure for multivariate polynomials that is being added to Maple for a future release.
9xyz - 4yz - 6xyz - 8x - 5
You can see the new structure above. It has the variables (x,y,z) followed by monomials and coefficients, where each monomial is encoded as a machine word. The encoding includes the total degree, so for xyz we encode the values (a+b+c,a,b,c), which is (a+b+c) 2 + a 2 + b 2 + c 2 on a 32-bit computer. The terms are then sorted in the encoding, which is effectively graded lexicographical order.
The purpose of this structure is to make computations with multivariate polynomials "as fast as possible". This is joint work underway with Dr. Michael Monagan at Simon Fraser University and with Maplesoft. Consider the following common polynomial operations:
All of these operations previously had to walk the entire Maple dag structure, but with this structure they are effectively constant time. Note that the number of variables is both limited and small. Further improvements are expected for other commands. For example, coeffs(f,x), diff(f,x), and eval(f,x=number) can all be computed in one pass through f without sorting the result. That is not possible in the Maple dag structure.
We expect the new structure to greatly improve Maple's ability to cope with very large polynomials, however fundamental changes like this take time. Maple 13 already uses this structure to internally to multiply and divide polynomials mod p. You can see the effect by comparing a computation mod p to a computation over the integers.
# Maple 13
p := prevprime(kernelopts(maximmediate));
f := expand((1+x+y+z)^20) mod p:
g := f+1:
time(assign('h'=Expand(f*g) mod p)); # new code
time(Divide(h,f,'q') mod p); # new code
time(assign('h'=expand(f*g))); # Maple dag
time(divide(h,f,'q')); # Maple dag
In Maple 14 this is extended to the integer case (expand and divide), which provides a big boost to Maple's overall polynomial performance. This still involves conversions to and from Maple's dag structure however. The full benefits will only be seen when we make this new structure the default.