64 Bit Monomial-Multiply Versions of Fermat

November 9, 2013

64 bit versions of Fermat. copyright Robert H. Lewis, 1995 - 2013.


Version 5.1, November 9, 2013. Version 5.0, October 25, 2013.

Introduction

Version 5.0 of 64 bit Fermat implements a new method for multiplication of multivariate polynomials. This provides very impressive speedups in real problems, often 40 - 70%, even more in some test cases. The method can be modified by the user command &_o.


Version 5.1 has a few bug fixes and small improvements, mostly regarding when to use the old standard multiplication.


Detailed Description

Though this way of multiplying polynomials actually has a long history, I am indebted to Michael Monagan of Simon Fraser University for bringing it to my attention at the 2013 ECCAD conference in Maryland. He has described the idea in a preprint available here.


Basically, the idea is to store each term, or monomial, of a multivariate polynomial in a single node instead of storing the polynomial as a recursive structure of nodes at one level pointing to nodes at a lower level.


For example, consider (x^2 + y^2 + z^3 +1)^3. This polynomial has 20 terms. As a recursive list structure it is
z^9 + (3y^2 + 3x^2 + 3)z^6 + (3y^4 + (6x^2 + 6)y^2 + 3x^4 + 6x^2 + 3)z^3 + y^6 + (3x^2 + 3)y^4 + (3x^4 + 6x^2 + 3)y^2 + x^6 + 3x^4 + 3x^2 + 1
This means that at the highest level, z, there are four nodes, one each for z^9, z^6, z^3, and z^0. Each node has a pointer to its coefficient, a polynomial at the next level down, for y. The coefficient of z^6 is the polynomial 3y^2 + 3x^2 + 3. That has two nodes, one for y^2 and for one for y^0. The coefficient of y^0 is the polynomial 3x^2 + 3, which is one more level down, and so forth. (Eventually there is a pointer for the numerical coefficient, stored in a different kind of node.) This has always been the storage structure of polynomials in Fermat. It has many advantages, one of which is that any exponent may be easily accommodated. Another is that multiplication is easily written as a recursive program; each level is handled in the same way as all others.


However, the disadvantage of the recursive structure is that there are many links to follow, and much space is devoted to storing these links.


In contrast, the monomial structure stores the above example as a list of twenty nodes, each of which contains all the exponents for that term. The polynomial is thought of as
z^9 + 3*y^2*z^6 + 3*x^2*z^6 + 3*z^6 + 3*y^4*z^3 + 6*x^2*y^2*z^3 + 6*y^2*z^3 + 3*x^4*z^3 + 6*x^2*z^3 + 3*z^3 + y^6 + 3*x^2*y^4 + 3*y^4 + 3*x^4*y^2 + 6*x^2*y^2 + 3*y^2 + x^6 + 3*x^4 + 3*x^2 + 1
For example, 6*x^2*y^2*z^3 is a node containing a pointer to the numerical coefficient, then the three exponents 2 2 3. These are stored as fields within a single long integer. For example, if one knew that only three-variable polynomials would ever be considered, one could use a 32 bit integer to store all three, with the exponent for x (lowest precedence) going in spots 0-9, for y in spots 10-19, and for z in spots 20-29. The upper two spots 30 and 31 would be unused.


The advantage is that multiplication of terms is now very easy -- just add the nodes containing the three exponents, a single 32 bit integer addition! (Of course, one multiplies the coefficients too.) This is fine as long as one does not encounter an exponent as large as 2^10 = 1024. Then disastrous overflow would occur, perhaps crashing Fermat, or worse, introducing undetected nonsense.


In Fermat we expect to use far more than three polynomial variables. I have therefore implemented this idea with 128 bit long long integers instead of 32 bit integers. The latest gcc C compiler has builtin addition and subtraction of such long integers. The user may choose the "monolength", the number of bits allocated to each variable's exponent (it was 10 in the example above) with the command &_o. Fermat will then compute and display monovars = 128 div monolength, mononum = 2^monolength, and monogap = 128 - monolength*monovars. The default monolength is 9.


If after the &_o prompt you enter a value > 30 or < 5, the monomial multiplication method is disabled; multiplication will proceed recursively. Since legal values are between 5 and 30, possible monovars are between 25 and 4. It is probably unwise to use a monolength < 7, so the most polynomial variables that can reasonably be multiplied by this method is 18. Therefore, one may well have a situation where one has more polynomial variables than mononvars. This is fine: Fermat begins the standard, recursive method, and switches when the depth of recursion allows it (depth remaining <= monovars) for multiplying coefficients.


There is no guarantee against overflow, other than the wisdom of the user. If overflow is suspected, there will be warning messages displayed. If enough messages occur, an error occurs.


There is no change in addition or subtraction of polynomials. The user's variables are still stored in the recursive structure. When the need to multiply arises, they are converted to the new format, multiplication is done, and the answer converted to the recursive form for storage.


What is all this for? Speed! Unless one has only two variables, or one of the multiplicands has few terms, say < 60, a noticeable improvement in speed results. On Michael Monagan's test cases in the above preprint, I do not see quite the speedups he reports, but the improvement is often remarkable. More important, on real problems taking many minutes to complete involving thousands of multiplications of polynomials with six or more variables, I routinely see 30 - 80% improvements.


Be aware: Again, there is no guarantee against overflow, or that you will be warned of it. You may get a "Fermat error", especially "bad exps in fast_poly divide".


Summary

A new fast multiplication method for polynomials is implemented in Fermat 5.x for 64 bit Mac and Linux. Adjust it with &_o. Experiment.








THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
NO ONE MAY USE THIS SOFTWARE AS PART OF ANY BUSINESS OR FOR-PROFIT ENDEAVOR WITHOUT THE EXPRESS WRITTEN PERMISSION OF ROBERT H. LEWIS. ROBERT H. LEWIS RETAINS THE COPYRIGHT TO ANY PROGRAM THAT INCLUDES THIS SOFTWARE.

(Thanks WebCounter!)