Enumerate grid points on 2D plane with descending order of (x * y)

Question

Given N > 0 and M > 0, I want to enumerate all (x, y) pairs such that 1 <= x <= N and 1 <= y <= M in descending order of (x * y). An example: given N = 3 and M = 2, the enumeration sequence should be:

1. (3, 2) -- 3 * 2 = 6
2. (2, 2) -- 2 * 2 = 4
3. (3, 1) -- 3 * 1 = 3
4. (2, 1) -- 2 * 1 = 2
5. (1, 2) -- 1 * 2 = 2
6. (1, 1) -- 1 * 1 = 1

The order of (2, 1) and (1, 2) could be swapped. One obvious way is to list them all, insert into a vector<pair<int, int> >, and call std::sort() with my own comparison function. However, since N and M may be large, and most of the time I only need the first few terms of the sequence, I hope there could be some smarter way that generates such a sequence instead of generate them all-and-sort, which requires as many as N*M array elements.

Update: I forgot to mention that although most of the time I only need the first few terms, the number of terms required is unknown before enumeration.

Answer 1

If you're just looking to save on space while retaining the time as more or less equal, you can count on the fact that each successively smaller element must be adjacent (in the 2-D grid you alluded to) to one of the elements you've already encountered. (You can prove this with induction, it's not particularly difficult. I'm going to assume for the rest of this that M>=N.)

The basic algorithm looks something like:

Start with a list (Enumerated Points) containing just the maximum element, M*N
Create a max heap (Candidates) containing (M-1),N and M,(N-1).
Repeat this:
    1.Pick the largest value (x,y) in Candidates and append it to Enumerated Points
    2.Add (x-1),y and x,(y-1) to Candidates if they are not there already

You can repeat this as long as you want more elements in Enumerated Points. The max size of Candidates should be M+N so I think this is O(k log(M+N)) where k is the number of points you want.

ADDENDUM: The matter of avoiding duplicate is not entirely difficult but is worth mentioning. I will assume in this algo that you lay your grid out so that numbers go down as you move down and right. Anyway, it goes like this:

At the beginning of the algorithm create an array (Column Size) which has one element for each column. You should make this array contain the number of rows in each column which are part of the list of enumerated points.

After you add a new element and update this array, you will check the size of the column on either side in order to decide if the grid squares to the immediate right and below of this new enumerated point are already in the candidates list.

Check the size of the column to the left- if it's larger than this one, you don't need to add the element below your new enumerated point.

Check the size of the column to the right- if it's one less than the same size of this column, you don't need to update the element to the right of this one.

To make this obvious, let's look at this partially completed chart for M=4, N=2:

4  3  2  1
*  *  *  2 |2
*  3  2  1 |1

The elements (4,2), (3,2), (2,2) and (4,1) are already in the list. (The first coordinate is M, the second is N.) The Column Size array is [2 1 1 0] since this is the number of items in each column that are in the Enumerated Points list. We are about to add (3,1) to the new list- We can look to the column size to the right and conclude that adding (2,1) isn't needed because the size of the column for M=2 is larger than 1-1. The reasoning is pretty clear visually - we already added (2,1) when we added (2,2).

Answer 2

Here's an O(K logK) solution, where K is the number of terms you want to generate.
Edit: Q keeps only one copy of each element; an insert fails if the element is already present.

priority_queue Q
Q.insert( (N*M, (N,M)) ) // priority, (x,y) 
repeat K times:
    (v, (x,y)) = Q.extract_max()
    print(v, (x,y))
    Q.insert( (x-1)*y, (x-1,y) )
    Q.insert( x*(y-1), (x,y-1) )

This works because before visiting (x,y) you must visit either (x+1,y) or (x,y+1). The complexity is O(KlogK) since Q at most 2K elements are pushed into it.

Answer 3

This is effectively equivalent to enumerating the prime numbers; the numbers you want are all the numbers that aren't prime (except for all those that have x or y equal to 1).

I'm not sure there's a method of enumerating primes that's going to be quicker than what you're already proposing (at least in terms of algorithmic complexity).

Answer 4

Because you mention that most of the time you need the first few terms of the sequence; after generating them all you don't need to sort them all to find these first few terms. You can use a Max Heap depending on number of terms that you want, say k. So if a heap is of size k (<< N && << M) then you can have largest k terms after nlogk which is better than nlogn for sorting.

Here n = N*M

Answer 5

A dummy approach that loops from NxM to 1 searching for pairs that when multiplied produce the current number:

#!/usr/bin/perl

my $n = 5;
my $m = 4;

for (my $p = $n * $m; $p > 0; $p--) {
    my $min_x = int(($p + $m - 1) / $m);
    for my $x ($min_x..$n) {
        if ($p % $x == 0) {
            my $y = $p / $x;
            print("x: $x, y: $y, p: $p\n");
        }
    }
}

For N=M, complexity is O(N³) but memory usage is O(1).

Update: Note that the complexity is not as bad as it seems because the number of elements to generate is already N². For comparison, the generate-all-the-pairs-and-sort approach is O(N²logN) with O(N²) memory usage.

Answer 6

Here's an algorithm for you. I'll try to give you an English description.

In the rectangle we're working with, we can always assume that a point P(x, y) has a greater "area" than the point below it P(x, y-1). So when looking for points of maximum area, we only need to compare the topmost untaken point in each column (i.e. for each possible x). For example, when considering the pristine 3 x 5 grid

5 a b c
4 d e f
3 g h i
2 j k l
1 m n o
  1 2 3

we really only need to compare a, b, and c. All the other points are guaranteed to have less area than at least one of those. So build a max heap that only contains the topmost point in each column. When you pop a point out of the heap, push in the point that's directly below it (if that point exists). Repeat until the heap is empty. This gives you the following algorithm (tested, but it's in Ruby):

def enumerate(n, m)
    heap = MaxHeap.new
    n.times {|i| heap.push(Point.new(i+1, m))}

    until(heap.empty?)
        max = heap.pop
        puts "#{max} : #{max.area}"

        if(max.y > 1)
            max.y -= 1
            heap.push(max)
        end
    end
end

This gives you a running time of O((2k + N) log N). Heap operations cost log N; we do N of them when building the initial heap, and then 2k when we pull out the k points of maximum area (2 assuming each pop is followed by the push of the point below it).

It has the additional advantage of not having to build all your points, unlike the original proposal of building the entire set and then sorting. You only build as many points as necessary to keep your heap accurate.

And finally: Improvements can be made! I haven't done these, but you could get even better performance with the following tweaks:

1. Only descend to y = x in each column instead of y = 1. To generate the points that you're no longer checking, use the fact that the area of P(x, y) is equal to the area of P(y, x). Note: If you use this method, you'll need two versions of the algorithm. Columns work when M >= N, but if M < N you'll need to do this by rows instead.
2. Only consider the columns that could possibly contain the maximum. In the example I gave, there's no reason to include a in the heap from the beginning, since it's guaranteed to be less than b. So you only need to add columns to the heap when the top points of their neighbour columns get popped.

And that turned into a small essay... Anyways - hope it helps!

Edit: The full algorithm, incorporating both improvements I mentioned above (but still in Ruby, 'cause I'm lazy). Note that there's no need for any extra structures to avoid inserting duplicates - unless it's a "top" point, each point will only insert another point in it's row/column when taken.

def enumerate(n, m, k)
    heap = MaxHeap.new
    heap.push(Point.new(n, m))
    result = []

    loop do
        max = heap.pop
        result << max
        return result if result.length == k

        result << Point.new(max.y, max.x) if max.x <= m && max.y <= n && max.x != max.y
        return result if result.length == k

        if(m < n) # One point per row
            heap.push(Point.new(max.x, max.y - 1)) if max.x == n && max.y > 1
            heap.push(Point.new(max.x - 1, max.y)) if max.x > max.y
        else # One point per column
            heap.push(Point.new(max.x - 1, max.y)) if max.y == m && max.x > 1
            heap.push(Point.new(max.x, max.y - 1)) if max.y > max.x
        end
    end
end

Answer 7

I got it!

Consider the grid as a set of M columns where every column is a stack containing the elements from 1 at the bottom to N at the top. Every column is tagged with its x coordinate.

The elements inside every column stack are ordered by its y value and so also by x*y as x has the same value for all of them.

So, you just have to go picking the stack that has the bigger x*y value at its top, pop it and repeat.

In practice you will not need stacks, just the index of the top value and you can use a priority queue to get the column with the bigger x*y value. Then, decrement the value of the index and if it is bigger than 0 (indicating that the stack has not been exhausted) reinsert the stack on the queue with its new priority x*y.

The complexity of this algorithm for N=M is O(N²logN) and its memory usage O(N).

Update: Implemented in Perl...

use Heap::Simple;

my ($m, $n) = @ARGV;

my $h = Heap::Simple->new(order => '>', elements => [Hash => 'xy']);
# The elements in the heap are hashes and the priority is in the slot 'xy':

for my $x (1..$m) {
    $h->insert({ x => $x, y => $n, xy => $x * $n });
}

while (defined (my $col = $h->extract_first)) {
    print "x: $col->{x}, y: $col->{y}, xy: $col->{xy}\n";
    if (--$col->{y}) {
        $col->{xy} = $col->{x} * $col->{y};
        $h->insert($col);
    }
}

Answer 8

In Haskell it produces the output immediately. Here's an illustration:

         -------
        -*------
       -**------
      -***------
     -****------
    -*****------
   -******------
  -*******------

Each starred point produces both (x,y) and (y,x). The algorithm "eats into" this thing from the top right corner, comparing the top elements in each column. The length of the frontier is never more than N (we assume N >= M).

enuNM n m | n<m = enuNM m n                    -- make sure N >= M
enuNM n m = let
    b = [ [ (x*y,(x,y)) | y<- [m,m-1..1]] | x<-[n,n-1..m+1]]
    a = [ (x*x,(x,x)) : 
          concat [ [(z,(x,y)),(z,(y,x))]       -- two symmetrical pairs,
                           | y<- [x-1,x-2..1]  --  below and above the diagonal
                           , let z=x*y  ] | x<-[m,m-1..1]]
 in
    foldi (\(x:xs) ys-> x : merge xs ys) [] (b ++ a)

merge a@(x:xs) b@(y:ys) = if (fst y) > (fst x) 
                            then  y : merge  a ys 
                            else  x : merge  xs b
merge a [] = a 
merge [] b = b

foldi f z []     = z
foldi f z (x:xs) = f x (foldi f z (pairs f xs))

pairs f (x:y:t)  = f x y : pairs f t
pairs f t        = t

foldi builds a skewed infinitely deepening tree serving as a heap, joining all the producer streams, each for each x, which are created already sorted in descending order. Since all the initial values of producer streams are guaranteed to be in decreasing order, each initial value can be popped without comparison, allowing the tree to be built lazily.

The code for a produces the pairs above diagonal line using the corresponding pairs from below the diagonal line (under the assumption N >= M, for each (x,y) where x <= M & y < x, (y,x) is also to be produced.)

It should be practically O(1) for each of the few first values produced which are very near the top of the tree of comparisons.

Prelude Main> take 10 $ map snd $ enuNM (2000) (3000)
[(3000,2000),(2999,2000),(3000,1999),(2998,2000),(2999,1999),(3000,1998),(2997,2
000),(2998,1999),(2999,1998),(2996,2000)]
(0.01 secs, 1045144 bytes)

Prelude Main> let xs=take 10 $ map (log.fromIntegral.fst) $ enuNM (2000) (3000)
Prelude Main> zipWith (>=) xs (tail xs)
[True,True,True,True,True,True,True,True,True]

Prelude Main> take 10 $ map snd $ enuNM (2*10^8) (3*10^8)
[(300000000,200000000),(299999999,200000000),(300000000,199999999),(299999998,20
0000000),(299999999,199999999),(300000000,199999998),(299999997,200000000),(2999
99998,199999999),(299999999,199999998),(299999996,200000000)]
(0.01 secs, 2094232 bytes)

We can assess the empirical run-time complexity:

Prelude Main> take 10 $ drop 50000 $ map (log.fromIntegral.fst) $ enuNM (2*10^8)
 (3*10^8)
[38.633119670465554,38.633119670465554,38.63311967046555,38.63311967046554,38.63
311967046554,38.63311967046553,38.63311967046553,38.633119670465526,38.633119670
465526,38.63311967046552]
(0.17 secs, 35425848 bytes)

Prelude Main> take 10 $ drop 100000 $ map (log.fromIntegral.fst) $ enuNM (2*10^8
) (3*10^8)
[38.63311913546512,38.633119135465115,38.633119135465115,38.63311913546511,38.63
311913546511,38.6331191354651,38.6331191354651,38.633119135465094,38.63311913546
5094,38.63311913546509]
(0.36 secs, 71346352 bytes)

*Main> let x=it
*Main> zipWith (>=) x (tail x)
[True,True,True,True,True,True,True,True,True]

Prelude Main> logBase 2 (0.36/0.17)
1.082462160191973     -- O(n^1.08) for n=100000 values produced

This can be translated into e.g. Python using generators for Haskell streams in a strightforward manner as can be seen here.