How to get the (x,y) coordinates given a index when x < y?

Question

In the past I've often used loops of the following kind (Haskell example):

  upperBoundToTuples :: Int -> [(Int, Int)]
  upperBoundToTuples n = [(x,y) | x <- [0..n], y <- [x+1..n]]

The above code produces tuples of range (0,1)..(n,n) where for all x < y.

I was wondering if it there was an efficient way of getting those (x,y) indices given a single index? Possible applications include optimization problems on the GPU where loops are not allowed and each thread only gets an index.

Also if it is possible for the 2D case, could such an algorithm be generalized to multiple dimensions?

Answer 1

You're asking for a bijection from [0, N(N+1)/2) to pairs (x, y) with 0 <= x < y <= N.

Here's one simple way to define it (in pseudocode, but should be trivial to convert to Haskell):

x0, y0 = i / (N + 1), i % (N + 1)
if x0 < y0 then result = (x0, y0)
else result = (N - 1 - x0, N - y0)

Here's a visualisation of the function for N=6. The map is laid out in a table with rows of length N+1=7, with the first row representing the value of the function for i=0 to 6, the next i=7 to 13 and so on. If you look very closely and carefully, you can see that things above the leading diagonal map to their own location in the table, and things on or below the diagonal map rotationally to the later entries.

5,6  0,1  0,2  0,3  0,4  0,5  0,6
4,6  4,5  1,2  1,3  1,4  1,5  1,6
3,6  3,5  3,4  2,3  2,4  2,5  2,6

And here's the opposite of this visualisation: a table T of size (N+1) by (N+1) with T[x, y] = i where i is mapped to (x, y) by the function above.

 -   1   2   3   4   5   6
 -   -   9  10  11  12  13
 -   -   -  17  18  19  20
 -   -   -   -  16  15  14
 -   -   -   -   -   8   7
 -   -   -   -   -   -   0
 -   -   -   -   -   -   -

Higher dimensions

This method can probably be made to work in higher dimensions, but I don't immediately see how. As an alternative, here's a simple but somewhat inefficient method that does work in arbitrary dimensions.

First, note there's choose(N + 1, k) increasing sequences of length k from the numbers from 0 to N (where choose(N, k) is the binomial coefficient). Of those, choose(N, k - 1) of them end with N. That gives this recursive function that generates the sequences in descending colexicographical order (again in pseudocode):

sequence(N, k, index)
    = []                                           if k == 0
    = sequence(N - 1, k - 1, index) + [N]          if index < choose(N, k - 1)
    = sequence(N - 1, k, index - choose(N, k - 1)) otherwise

Here's, sequence(5, 3, index) for index between 0 and 19:

 0 -> [3, 4, 5]
 1 -> [2, 4, 5]
 2 -> [1, 4, 5]
 3 -> [0, 4, 5]
 4 -> [2, 3, 5]
 5 -> [1, 3, 5]
 6 -> [0, 3, 5]
 7 -> [1, 2, 5]
 8 -> [0, 2, 5]
 9 -> [0, 1, 5]
10 -> [2, 3, 4]
11 -> [1, 3, 4]
12 -> [0, 3, 4]
13 -> [1, 2, 4]
14 -> [0, 2, 4]
15 -> [0, 1, 4]
16 -> [1, 2, 3]
17 -> [0, 2, 3]
18 -> [0, 1, 3]
19 -> [0, 1, 2]

Answer 2

We may equivalently consider [(x,y) | x<-[0..n], y<-[0..x-1]]. This list has length

ℓ_n = _x=0∑ⁿ x = ^n·(n+1)/₂.

Hence we can get, to a given ℓ, the nearest lower n through

2·ℓ_n = n·(n+1) = n² + n

n_~ = -½ ± √(¼ + 2·ℓ_n)

In particular, for a given index i,

n_i− = ⌊-½ ± √(¼ + 2·i)⌋

is the x-length of the last fully completed triangle. Thus, the index i lies in the row n_i−+1. That triangle had an area of

ℓ_ni− = ^{n_i−·(n_i−+1)}/₂

which we therefore need to subtract from i to get the remainder index (in y-direction). This gives rise to the definition

lowerTriangularTuple :: Int -> (Int,Int)
lowerTriangularTuple i = (nmin+1, i - (nmin*(nmin+1))`div`2)
 where nmin = floor $ -1/2 + sqrt(1/4 + 2 * fromIntegral i)

Example:

GHCi> lowerTriangularTuple <$> [0..30]
[(1,0),(2,0),(2,1),(3,0),(3,1),(3,2),(4,0),(4,1),(4,2),(4,3),(5,0),(5,1),(5,2),(5,3),(5,4),(6,0),(6,1),(6,2),(6,3),(6,4),(6,5),(7,0),(7,1),(7,2),(7,3),(7,4),(7,5),(7,6),(8,0),(8,1),(8,2)]