# De-linearize a 2d linearized grid and calculate distance between 2 points

In a simulation of mine, for which performance is paramount, I have the following bottleneck.

I have a huge 2D grid of dimensions width and height, whose dimensions are stored in linearized form, i.e. positions at x and y are stored as one integer x + y*width. Every iteration I have to get thousands of such linearized coordinates, de-linearize an calculate distance between 2d coordinates.

Currently, what I do is:

y1 = linPos1 / width; //int divided by int disregards remainders
x1 = linPos1 - width * y1;

y2 = linPos2 / width; //int divided by int disregards remainders
x2 = linPos2 - width * y2;

int xDiff = x2-x1;
int yDiff = y2-y1;
double distance = sqrt(xDiff*xDiff  + yDiff*yDiff);


How could I improve upon that? What is the most performant way to take two of such integer values, retrieve back x and y and then calculate the 2D distance between said points? I am open to all suggestions, including SIMD instructions.

• What happens if width=2, x=2, y=3? The value stored is 2+3*2. If you divide that by width, you'll get 4 as y, and then 8 - 2*4 = 0 as x. Isn't this a problem? – Andrew Apr 19 '17 at 15:15
• @Andrew I should have added that this is a square grid always. So that won't ever happen. – Louis15 Apr 19 '17 at 15:17
• So can we assume that $0 \leq x,y \leq \text{width}-1$? – Andrew Apr 19 '17 at 15:23
• @Andrew since the scale of the axis is indifferent for me, yes. I could easily adapt the code so the grid is bounded any way needed. – Louis15 Apr 19 '17 at 15:29
• One last question, do you need the euclidean metric, or would some other one do too? – Andrew Apr 19 '17 at 15:31

What you are doing is packing the two integer coordinates into a single integer, then unpacking them.

Let us assume that $\text{height} \leq \text{width}$, let w = width. Translate the lattice so that $0 \leq x,y\leq w-1$.

Using these assumptions, the packing can be done by writing the coordinates in base w. For example, if w = 10, x=4, and y=7, then the representation is w*y + x = 7*10 + 4 = 74.

Using our assumptions this is an injective function, i.e. to one pair $(x,y)$ there exists at most one $k$ such that $k=wy+x$, and each such $k$ uniquely determine the $(x,y)$. This means that the packing-unpacking is well-defined, which is great: same points will have the same representations, and a representation corresponds to exactly one point.

y = k / w;
x = k - w * y;


could be written a few different ways, e.g. using %,

y = k / w
x = k % w


or by using std::div.

Modern compilers will compile all of these to the same instruction: idivl, so you don't have to worry about this part, just don't forget to tell your compiler to optimize slightly. For example, g++ needs the -O1 flag.

So that's about as good as it gets, unless we can use some faster form of integer division(from now on, simply division).

If a number is in base $b$, then it's really easy to divide it by powers of $b$: just shift the numbers around; e.g. 1234/10^3 = 1, 1234 % 10^3 = 4.

That's all fine and well, but string operations aren't going to be faster than idivl. Luckily we don't need those, because by choosing $w=2^n$, we can leverage bit-shifts, and bitwise-and to do exactly the same for us, really fast.

Let's find the smallest $n \in \mathbb N$, such that $w \leq 2^n$. Let $w = 2^n$. Perhaps you won't use a few coordinates, but that's not really a problem. Now our division code will look like this.

y = k >> n
x = k & ((1 << n) - 1)


Here x will be the first nth bits of k, and y will be the rest. For example, if k = 1101, n=2, then
x = 1101 & ((1 << 2) - 1) = 1101 & (0100 - 0001) = 1101 & 0011 = 0001 = 01,
y = 1101 >> 2 = 0011 = 11.

Now let's see if the compilers are clever enough to do this for us automatically.

int func1(int a) {
int x = a / 1024;
int y = a % 1024;
return x+y;
}

int func2(int a) {
int x = a >> 10;
int y = a & ((1 << 10) - 1);
return x+y;
}


Let's see the assembly code generated from the code above, if we allow the compiler to optimize a bit.

_Z5func1i:
leal    1023(%rdi), %eax
testl   %edi, %edi
movl    %edi, %edx
cmovns  %edi, %eax
sarl    $31, %edx shrl$22, %edx
sarl    $10, %eax addl %edx, %edi andl$1023, %edi
subl    %edx, %edi
ret

_Z5func2j:
movl    %edi, %eax
andl    $1023, %edi shrl$10, %eax
ret


As you can see, besides some sign-checking, the two codes are identical. In fact, after specifying that we are working with unsigned integers, we get the same output:

    movl    %edi, %eax
andl    $1023, %edi shrl$10, %eax
ret


This means that we can keep our more readable operators: %, and /. Hurray!

Now let's look at the last part of your code, and let's change it to a more general formula.

double distance = metric(xDiff,yDiff);


Where metric is a function calculating a distance. In your case, metric corresponds to the Euclidean metric,

double metric(int dx, int dy) {
return std::sqrt(dx*dx+dy*dy);
}


Depending on your needs, you could use other metrics; as you are working on an integer lattice, I'd suggest you take a look at the taxicab metric. Its calculation is faster, and there is less room for numerical errors.

double metric(int dx, int dy) {
return std::abs(dx) + std::abs(dy);
}


Implementing these changes result in code, which is a constant time faster, than your original one.

While it was fun figuring out, I doubt this is the kind of speed up you were hoping for. I think optimizing other parts of your code could have greater benefits. Perhaps consider posting them as a separate question.

• Many thanks. I see what you did, the trick is actually in how de-linearize the 2D array. I like that, due to the flexibility. I will give it a try in a couple of hours and report back (I can't upvote due to lack of rep, otherwise I would) – Louis15 Apr 19 '17 at 16:26
• @Louis15 No problem. You can still accept the answer, if -after a while- this one is the best in your opinion. By the way, I've just now edited in additional explanation. – Andrew Apr 20 '17 at 13:17

You can obtain both the quotient y and the remainder x using just one call to std::div(linPos, width).

If width happens to be a power of 2, then you can quickly get the quotient using a right shift (>>), and the remainder using bitwise AND (&).

• I just want to add that g++ even only with the -O1 flag compiles OP's corresponding code, and std::div to the exact same instruction, namely idivl. – Andrew Apr 19 '17 at 18:09