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.

Answer 3

As far as my understanding goes, you cannot even hope for a solution with complexity O(N*M) because of the following logic.

Say the arrays are (a1, a2, a3, a4) and (b1, b2, b3, b4, b5)

The possible pairs shall be as follows:

a1b1 a1b2 a1b3 a1b4 a1b5
     a2b1 a2b2 a2b3 a2b4 a2b5
          a3b1 a3b2 a3b3 a3b4 a3b5
               a4b1 a4b2 a4b3 a4b4 a4b5

Now for each row, the pairs at the left shall always be picked up before the pairs on the right.

1. So 1st candidate is a1b1.
2. Then once a1b1 gets picked, the next candidates are a1b2 and a2b1.
3. Say a1b2 gets picked. Then the candidates are a1b3 and a2b1.
4. Now say a2b1 gets picked. Then the candidates are a1b3, a2b2 and a3b1.

Hence we see that as we move forward, the number of candidates for the position increase linearly.

Hence in the worst case, you will have to carry out comparisons of the order of N*M*M.

An O(N*Mlog(M*N)) approach. (where N = a.size() and M = b.size())

for ( int i = 0; i < a.size(); i++ )
  for ( int j = 0; j < b.size(); j++ )
    sumVector.push_back( a[i] + b[j] );

sort( sumVector ); // using merge sort or quicksort.

Answer 4

I'm hoping it's possible storing only a constant amount of state and only iterating over each list once.

I think this is impossible. I can't provide a rigorous mathematical proof, but consider this: your problem has two parts:

1) Generate all pairs from A and B. 2) Make sure they are in order of increasing sum.

Let's drop the second part to make it easier. The simplest possible way to implement this is

foreach a in A:
  foreach b in B:
    emit (a, b)

I hope you will agree that there can be no more efficient way of doing this. But here we already iterate over B length(A) times.

So the minimum time complexity for this is O(A*B), but you want O(A+B).

Answer 5

It's quite simple using an algorithm to find a pair of number in a sorted array that sum to a specific constant K, as described in this question.

Summary

The final solution would be of time complexity O((M+N)*(M+N)), which is at most 4 times worse than the lowerbound of M*N (just outputting all the pairs). So constant factor is only 4.

Edit: Oops, it's not O((M+N)*(M+N)) as I thought, apparently, it's O(K*(M+N)), where K is the highest sum in the two arrays. I'm not sure whether this algorithm can be improved, but it seems that that solution will be similar to the Fast Fourier Transform method described by Peter de Rivaz.

Sum in sorted array algorithm

In that algorithm, we set a lower pointer to 0, and higher pointer to the end of the array. Then if the sum at these two positions is larger than K, we decrease higher pointer. If it's lower, we increase the lower pointer. This works because at any iteration, arr[low] is the lowest element so far from which an answer might be derived. Similarly, arr[high] is the highest. So if we take the lowest and highest elements, and the sum is larger than K, we know that any other combination with arr[high] will be larger than K. So arr[high] can't be part of any solution. Therefore we can remove it from our array (this is achieved by reducing the high pointer).

Applying that to your problem

Applying this for your problem, the idea is, we iterate over the possible sum, that is from 0 to A[len(A)-1]+B[len(B)-1]. And for each possible sum, we run the above algorithm. For your problem, we set lower pointer to array A, and higher pointer to array B.

For the original algorithm, it will break as soon as it finds a pair which sum to the constant. For your problem, you will increase ptr_A and decrease ptr_B each by 1. This works because your array are strictly increasing. So if we find that A[ptr_A]+B[ptr_B]==K, then A[ptr_A]+B[low_B]<K for all low_B<ptr_B and A[high_A]+B[ptr_B]>K for all high_A>ptr_A.

And this will find all pairs, because for each possible sum K, it will find all pairs that sum to K, and we iterate over all possible sum K.

As a bonus, this algorithm will sort the output based on increasing value in list A (you can sort based on increasing value in list B by interchanging the pointers), and that we don't need random access to the array.

Code in Python

Implementation in Python:

def pair_sum(A,B):
    result = []
    max_sum = A[-1]+B[-1]
    for cur_sum in range(max_sum+1): # Iterate over all possible sum
        ptr_A = 0        # Lower pointer
        ptr_B = len(B)-1 # Higher pointer
        while True:
            if A[ptr_A]+B[ptr_B]>cur_sum:
                ptr_B -= 1
            elif A[ptr_A]+B[ptr_B]<cur_sum:
                ptr_A += 1
            else:
                result.append((A[ptr_A],B[ptr_B]))
                ptr_A += 1
                ptr_B -= 1
            if ptr_A==len(A):
                break
            if ptr_B==-1:
                break
    return result

def main():
    print pair_sum([0,1,2],[1,4])
    print pair_sum([0,1],[0,1])
    print pair_sum([0,1,3],[1,2])

if __name__=='__main__':
    main()

will print:

[(0, 1), (1, 1), (2, 1), (0, 4), (1, 4), (2, 4)]
[(0, 0), (0, 1), (1, 0), (1, 1)]
[(0, 1), (0, 2), (1, 1), (1, 2), (3, 1), (3, 2)]

as desired.

Answer 6

So let's take two 'scares arrays' sA and sB, that contains only non null degree/coefficients figures as described in the original post.
example :

A   =  x^5 + 3*x^2 + 4
sA  = [ 5, 1, 2, 3, 0, 4 ]

B = 2*x^6 + 5*x^3 + 8*x
sB = [ 6, 2, 3, 5, 1, 8]

What i suggest is to do the operation the way we would do it by hand, so they take a time in m*n, where m,n are the number of non-null coefficients, and not in p*q, where p and q are the degrees of A, B.
Since you said m and n are small, m*n is no big deal.
To store the coefficients while computing, use either a sparse array (might be costly) or a hash table. the index or key is the degree, while the value is the corresponding coefficient.
Here an example of implementation in javascript, using a hash table :

http://jsbin.com/ITIgokiJ/2/edit?js,console

one example :

A     =   "x^5+3x^2+4"
B     =   "2x^6+5x^3+8x^1"
A * B =   "2x^11+11x^8+16x^6+15x^5+44x^3+32x^1"

the code :

function productEncodedPolynoms( sA, sB) {
   var aIndex = 0 ;
   var bIndex = 0 ;
   var resIndex = 0 ;
   var resHash = {} ;
   // for loop within sA, moving 2 items at a time
   for (aIndex = 0; aIndex < sA.length ; aIndex+=2) {
       // for loop within sB, moving 2 items at a time
       for (bIndex = 0; bIndex < sB.length ; bIndex+=2 ) {
           resIndex = sA[aIndex]+sB[bIndex] ;
           // create key/value pair if none created
           if (resHash[resIndex]===undefined)  resHash[resIndex]=0;
           // add this product to right coefficient
           resHash[resIndex] += sA[aIndex+1]*sB[bIndex+1];
       }
   } 
   // now unpack the hash into an encoded sparse array
   // get hash keys
   var coeff = Object.keys(resHash);
   // sort keys in reverse order
   coeff.sort(reverseSort);
   encodedResult = [];
   for (var i=0; i<coeff.length; i++ ) {
     if (resHash[coeff[i]]) {
       encodedResult.push(+coeff[i]);          // (+ converts to int)
       encodedResult.push(+resHash[coeff[i]]);
     }
   }
   return encodedResult;
}

example :

sA = [ 5, 1, 2, 3, 0, 4 ] ;

sB = [ 6, 2, 3, 5, 1, 8] ;

sAB = productEncodedPolynoms ( sA, sB );

utility to print result :

function printEncodedArray(sA) {
  res='';
  for (var i=0; i<sA.length; i+=2) {
    if (sA[i+1]) {
        if (sA[i+1] != 1 || sA[i]==0) res+=sA[i+1];
        if (sA[i]!=0) res+='x^'+sA[i];
        if (i!=sA.length-2) res+='+';
    }
  }
  return res;
}

// utilities
function reverseSort(a,b) { return b-a ; }

Answer 7

Why wouldn't you just generate the unsorted 2-d array then use a quicksort on it?

ED: Looks like if you generated the final array marginally intelligently (using a comparative algorithm during initial iteration and generation of the array) you could later use a more specific sorting algorithm (e.g. smoothsort) and end up with performance approaching O(2(m*n)) with a worst-case of O(m*n + (m*n) log (m*n))

ED2: I'm absolutely certain you could write an algorithm to do this in O(m*n). What you want to do is generate the array cleanly. I don't have time to write the pseudocode, but if you were to 1. Set up a minimum-iterator, maximum-iterator, and current-iterator for each of the arrays, as well as a current-max-sum for them. 2. Generate the first pair 3. Iterate one variable, compare it to the other possible third pair, (saving that in a temp variable) and either save it to the output or save the other array to the output and start again.

The way I visualize it is as two rows of blocks with three colors: red for completely finished, blue for 'highest so far' and unfinished. You work on one of the rows of blocks at a time, generating the next value for the result array, always comparing for the combinatorial of your baseline on the current side being lower than the currently-processing sum. Each time the sum is lower, you set the new lowest result to the result you just calculated, add the previously lowest result to your (already sorted) output array, and start adding up using the previous low (and its row) as a baseline. Never going back into the red blocks, and swiching which side is blue each time you have found a new opposite-side low.

You should never be calculating a single result twice, and come out with a sorted array.

It's essentially a floating-top bubble sort. You know the calculated value is higher than your current sum, and until it's not you process lower registers in the array. When it's no longer higher, you switch your top value up and go back to the point where you were iterating through the array whose members are currently smaller. EX:

A: (1 5 10)

B: (1 6 9)

(1, 1) -> (1, 6) is higher than (1, 5) so add (1, 5) make (1, 6) highest

(5, 6) is higher than (1, 6) so add (1, 6) make (5, 6) highest move up to (1, 9)

(1, 9) lower than (5, 6) so [add (1, 9) and increase finished on A]

Current array : (1, 1), (1, 5), (1, 6), (1,9)

(1, 10) equal to (5,6) so add both, increase finished on B, start with (5,9)

(5, 9) higher than (10, 1) so add (10, 1) make (5, 9) highest move up to (10, 6)

(10, 6) higher than (5, 9) so....

Anyways. If you set it up right, you have a single comparison made per addition to the final array, and get to build it on the fly without branching. If you don't need the final array in memory anywhere, the FFT's might go faster, but this should work.