3
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ 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? \$\endgroup\$ – Andrew Apr 19 '17 at 15:15
  • 1
    \$\begingroup\$ @Andrew I should have added that this is a square grid always. So that won't ever happen. \$\endgroup\$ – Louis15 Apr 19 '17 at 15:17
  • \$\begingroup\$ So can we assume that \$ 0 \leq x,y \leq \text{width}-1\$? \$\endgroup\$ – Andrew Apr 19 '17 at 15:23
  • \$\begingroup\$ @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. \$\endgroup\$ – Louis15 Apr 19 '17 at 15:29
  • \$\begingroup\$ One last question, do you need the euclidean metric, or would some other one do too? \$\endgroup\$ – Andrew Apr 19 '17 at 15:31
3
\$\begingroup\$

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.

Your code,

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
    addl    %edi, %eax
    ret

_Z5func2j:
    movl    %edi, %eax
    andl    $1023, %edi
    shrl    $10, %eax
    addl    %edi, %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
    addl    %edi, %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.

\$\endgroup\$
  • \$\begingroup\$ 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) \$\endgroup\$ – Louis15 Apr 19 '17 at 16:26
  • \$\begingroup\$ @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. \$\endgroup\$ – Andrew Apr 20 '17 at 13:17
0
\$\begingroup\$

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 (&).

\$\endgroup\$
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – Andrew Apr 19 '17 at 18:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.