Sparse polynomial division using a heap

Abstract

In 1974, Johnson showed how to multiply and divide sparse polynomials using a binary heap. This paper introduces a new algorithm that uses a heap to divide with the same complexity as multiplication. It is a fraction-free method that also reduces the number of integer operations for divisions of polynomials with integer coefficients over the rationals. Heap-based algorithms use very little memory and do not generate garbage. They can run in the cpu cache and achieve high performance. We compare our C implementation of sparse polynomial multiplication and division with integer coefficients to the routines of the Magma, Maple, Pari, Singular and Trip computer algebra systems.

Parallel sparse polynomial division using heaps

“…This is not optimal when #q < #g. In [13] we present a sequential division algorithm that does O(#f + #q#g log min(#q, #g)) comparisons. However our divisor heap is run on subproblems with #g/p by #q terms where p is the number of threads, so the threshold becomes easier to meet as the number of threads increases.…”

“…The strip height s is derived from the number of terms in g, refer to Section 2.3 for details. The threads merge terms from left to right in the style of a divisor heap of Monagan and Pearce [13]. Each iteration of the main loop extracts all of the products gi · qj with the largest monomial, multiplies their coefficients to compute a sum of like terms, and inserts their successors gi · qj+1 into the heap to set up the next iteration of the algorithm.…”

“…We can safely drop the products missing qj+1 as long as they are reinserted before they could be merged. For example, in our algorithm in [13] we set bits to indicate which gi have a product in the heap. When a new quotient term qj+1 is computed we check if g2 has a product in the heap and insert g2 · qj+1 if it does not, and when we insert gi · qj with i < #g, we also insert the next product for gi+1 if it is not already in the heap.…”

“…Given polynomials f and g with #f and #g terms, we construct f × g = P #f i=1 P #g j=1 fi · gj by creating, sorting, and merging all the products in parallel, entirely in the cache. We based the algorithm on Johnson's method [7] which we found to be a fast sequential approach in [12,13].…”

“…Then in Section 3 we present benchmarks of our implementation. We compare its performance and speedup to the sequential routine of [13], the parallel multiplication codes of [11], and the division routines of other computer algebra systems.…”