Iterating over pairs of numbers in an increasing order

Question

I'll explain where the question comes from at the bottom, but here's the statement. Suppose I have two lists of non-negative integers, which I'll write (A[0] ... A[n]) and (B[0] ... B[m]). They are strictly increasing, so A[i+1] > A[i] for all i and similarly for B. I want to collect all n * m pairs of elements in increasing order of their sum.

So, for example if A = (0 1 2) and B = (1 4), then I want to end up collecting ((0 1) (1 1) (2 1) (0 4) (1 4) (2 4)). If there is a tie, I don't care which order I collect the two elements in. For example, if A = (0 1) and B = (0 1), then I don't mind which of the mixed terms, (0 1) or (1 0), I pick up first.

Obviously, I'd like this to be reasonably efficient. I hope that it's possible in time asymptotic to m * n. Specifically, I'm hoping that the ordered input makes this problem strictly easier than the equivalent problem if I didn't know anything about the inputs. What I was thinking about when I first asked the question was the amount of state we have to store. I was hoping it was possible with a constant amount, but maybe this is unrealistic. (The things I've tried since have all failed!)

The code will actually be written in Lisp, but I think that the problem statement is pretty much independent of that. The inputs would most naturally come as singly-linked lists, but I will have to reverse them in advance anyway, so if random-access is relevant, I can make them arrays. In case it's relevant, I expect this to be mostly called on quite small lists, so a massive constant term / constant factor in runtime probably precludes a solution. (Although I'd be fascinated to hear about the algorithm idea!)

The background: I've been looking at the source code for Maxima, a computer algebra system and in particular at its code for the multiplication of two polynomials. Polynomials are represented in a "sparse format", so x^5 + x^2 + 2 might appear as (5 1 2 1 0 2), with descending exponents followed by their respective coefficients. To compute the product efficiently, what I really want to do is to collect together the degree zero terms, then the degree 1 terms etc. The current code avoids solving this problem by making a half-hearted stab at it for efficiency, then doing a sort of generic polynomial addition to deal with coefficients in an order it doesn't expect. I feel like we should be able to do better!

Answer 1

This problem differs only superficially from Sorting X + Y, a major irritant to computational geometers over the years. (See Joseph O’Rourke’s (open) Problem 41.) To summarize that link for implementers, the fastest known algorithm when the exponents can be added and compared only is O(m n log (m n)), and it’s the obvious one. If the exponents are bounded integers, then Peter’s Fourier approach applies. A lot of smart people have thought about this problem for quite a while, so I wouldn’t expect a better algorithm any time soon.

Answer 2

I wonder how sparse your polynomials are expected to be?

One option that is may be worth considering for the multiplication of dense polynomials is to compute the Fourier transform of the polynomials and multiply their Fourier coefficients together.

This allows you to multiply polynomials in O(nlogn) where n is the degree of the polynomial.

This would not be appropriate for a sparse polynomial such as (1+x^1000000)*(1+x^1000000) as n would be around 1000000 while the current algorithm would only need a few cycles.

A good explanation of this Fourier approach can be found in these lecture notes.