# Distance between two n-dimensional points (NASM)

I just finished writing a function that computes the distance between two n-dimensional points.

The original one was written in C and it's basically a translation of this formula:

$\text{dist}(x,y)=\sqrt{\sum_{i^1}^d{(x_i-y_i)^2}}$

Here's the C code for completeness sake:

void distanceC(float* source, float* destination, int dimensions, float* result){
// sqrt(sum((source-destination)^2, dimensions)))
int i;
for (i = 0; i < dimensions; i++){
*result += powf(source[i] - destination[i], 2);
}
*result = sqrtf(*result);
}


Finally, here's the NASM version that I wrote (I stripped some useless parts):

section .data

; masks used to set exceeding items to 0 when vectors are not multiple of 4

align 16
mask1:  dd      0xffffffff, 0x00, 0x00, 0x00
align 16
mask2:  dd      0xffffffff, 0xffffffff, 0x00, 0x00
align 16
mask3:  dd      0xffffffff, 0xffffffff, 0xffffffff, 0x00

section .text

Source      equ     8
Destination equ     12
Dimensions  equ     16
Result      equ     20

; ------------------------------------------------------------
; Distance between two points in n-dimensional space
; ------------------------------------------------------------
distance:
; ------------------------------------------------------------
; Entering the function
; ------------------------------------------------------------
push    ebp
mov     ebp, esp
sub     esp, 4
push    ebx
push    esi
push    edi

; ------------------------------------------------------------
; ------------------------------------------------------------
mov     esi, [ebp+Source]
mov     edi, [ebp+Destination]
mov     ecx, [ebp+Dimensions]
mov     eax, [ebp+Result]

; let's compute how many iterations do I need by doind Dimensions/4
; the result will tell how many packed iterations are needed, while the remainder will tell
; if I need we have exceeding items that require a final masked operation
push    eax
mov     edx, 0
mov     eax, ecx
mov     ebx, 4
div     ebx
mov     ecx, eax                ; ecx becomes the counter for the packed iterations
pop     eax

; now ecx contains the packed counter while edx contains the number of the exceeding items
; xmm0 = result
; xmm1 = source's components
; xmm2 = destination's components
; ebx  = both vector's offset
mov     ebx, 0
xorps   xmm0, xmm0              ; setting xmm0 to 0 (may not be needed)

; we may have a situation in which vector's size is less than 4, and this means
; that we don't have any packed iteration to do but we'll jump to the final
cmp     ecx, 0
je      .dU1

.loopP: movaps  xmm1, [esi+ebx*4]       ; let's get Source's next 4 items
movaps  xmm2, [edi+ebx*4]       ; and Destination's next 4 ones
subps   xmm1, xmm2              ; (Source-Destination)
mulps   xmm1, xmm1              ; (Source-Destination)^2
dec     ecx                     ; let's decrease ecx since we just did one iteration
haddps  xmm1, xmm1              ; let's sum the partial results that we have in xmm1
addps   xmm0, xmm1              ; let's add the new computed results to the total
jz      .loopC                  ; if ecx is 0 we're done with the packed iterations
add     ebx, 1                  ; let's move the offset to next elements
jmp     .loopP                  ; next iteration

; the packed loop is done, let's see if we have some exceeding items that are needing the masked iteration
.loopC  cmp     edx, 0
je      .endD                   ; edx = 0 means no masked iterations,

; if we didn't jump than we have to deal with the masked iteration
; we need to decide which mask to apply and we decide it based on edx's value

.dU1:   cmp     edx, 1
jne     .dU2                    ; if edx is not 1 than it may be 2
; if edx is 1 then we use mask1
jmp     .dU                     ; now that I have the correct mask I can jump to the actual execution

.dU2:   cmp     edx, 2
jne     .dU3
jmp     .dU

.dU3:   movaps  xmm7, [mask3]           ; if I'm here than I can only use the third mask since the previous checks failed

; let's apply the mask and complete the computation
.dU:    movaps  xmm1, [esi+ebx*4]       ; let's get Source's next 4 items
movaps  xmm2, [edi+ebx*4]       ; and Destination's next 4 ones

; now I can use the same instructions of the packed loop but just once
subps   xmm1, xmm2              ; (Source-Destination)
mulps   xmm1, xmm1              ; (Source-Destination)^2
haddps  xmm1, xmm1              ; let's sum the partial results that we have in xmm1
addps   xmm0, xmm1              ; let's add the new computed results to the total

; finally we need to apply the square root and copy the result where eax is pointing
.endD:  sqrtps  xmm0, xmm0
movss   [eax], xmm0

; ------------------------------------------------------------
; Exiting the function
; ------------------------------------------------------------

pop edi
pop esi
pop ebx
mov esp, ebp
pop ebp
ret


(You can test it by calling it from a C program, just add global distance)

The question is: is there any way to make it faster or just write it better?

## Wrong increment

First, there is a bug in your main loop. This line:

    add     ebx, 1


should be:

    add     ebx, 4


because you are operating on 4 floats at a time. Right now, you are doing array elements [0..3] following by [1..4]. On my computer the program crashed because it had a problem doing an unaligned load.

## Simplify divide by 4

This part where you divide ecx by 4:

    push    eax
mov     edx, 0
mov     eax, ecx
mov     ebx, 4
div     ebx
mov     ecx, eax                ; ecx becomes the counter for the packed iterations
pop     eax


can be simplified to this:

    mov     edx, ecx
and     edx, 3
shr     ecx, 2


Explanation: shr is the shift right instruction. When you shift right by 2 it is the same thing as dividing by 4. To get the remainder of something divided by 4, you just need to and that number with 3. C equivalent code:

    remainder = length & 3;
length >>= 2;


## Remove one jump from main loop

.loopP: movaps  xmm1, [esi+ebx*4]
movaps  xmm2, [edi+ebx*4]
subps   xmm1, xmm2
mulps   xmm1, xmm1
dec     ecx
jz      .loopC
jmp     .loopP

.loopC  cmp     edx, 0
je      .endD


Can be rewritten to remove one of the jumps at the end:

.loopP: movaps  xmm1, [esi+ebx*4]
movaps  xmm2, [edi+ebx*4]
subps   xmm1, xmm2
mulps   xmm1, xmm1
dec     ecx
jnz      .loopP

cmp     edx, 0
je      .endD


By the way, I also noticed another bug here. When you exit the loop you currently are skipping the add ebx, 4 line. But you need to advance ebx or else the remainder part below won't have the right value for ebx. This works correctly in my modified version because I always add to ebx inside the loop.

The part where you do a chain of compares to load one of three masks:

.dU1:   cmp     edx, 1
jne     .dU2
jmp     .dU

.dU2:   cmp     edx, 2
jne     .dU3
jmp     .dU



can be simplified to this:

.dU1:   shl     edx, 4
movaps  xmm7, [mask1 + edx - 16]


Explanation: Your masks are located in memory in consecutive order like an array of 128 bit values. So you can load them by index. Note that edx is 1..3 here. The shl instruction is shift left, and shifting left by 4 multiplies edx by 16. Then we load from mask1 + edx - 16. The minus 16 is because edx started at one instead of zero. Equivalent C code is something like this:

        uint8_t mask1[48] = { /* values */ };


## Stack push/pop

I'm not sure if this is part of your calling convention, but your prologue and epilogue have extraneous instructions:

    push    ebp
mov     ebp, esp
sub     esp, 4
push    ebx
push    esi
push    edi

pop     edi
pop     esi
pop     ebx
mov     esp, ebp
pop     ebp
ret


Could be changed to:

    push    ebp
mov     ebp, esp
push    ebx
push    esi
push    edi

pop     edi
pop     esi
pop     ebx
pop     ebp
ret


Furthermore, you are not utilizing eax. If you used eax instead of ebx everywhere in your function, you wouldn't have to save ebx because you wouldn't ever use it. You can load the result pointer into eax at the very end instead of at the beginning.

In general I liked all your comments. The one thing I thought you could add a comment for was that this double instruction:

        haddps  xmm1, xmm1


adds all 4 floats in xmm1 into the first float in xmm1. I wasn't familiar with the haddps instruction and I initially thought you made a typo (I thought a single haddps would do the whole job).

## Alignment required

I was able to use the function successfully, but when I modified my main(), the program crashed. The difference was that my modification caused the arrays I passed to distance() to be not aligned to 16 bytes. Since the movaps instruction expects 128 bit alignment, the arrays passed in need to have 128 bit alignment. I think you should switch to the movups instruction which handles unaligned floats.

• Thanks for your great feedback. I'm using add ebx, 1 because I've got movaps xmm1, [esi+ebx*4] so it's still moving each 4 items. I don't know which one is faster but I think that the code is correct. Removing the jump is a great idea, but can you please add some explanation to the division and the masking parts? I'd like to understand what I write. Finally, I'm using movaps because I'm aligning my data in the caller, so I'm sure that it won't lead to a segfault. – StepTNT May 14 '15 at 9:27
• @StepTNT A float is 4 bytes so ebx*4 will advance one float per one ebx. You need to advance 16 bytes per iteration, so add 4 to ebx. I'll try to add more explanation to the other parts. – JS1 May 14 '15 at 11:02
• I know that a float is 4 bytes, but I'm actually advancing on pointers, that's why ebx*4 moves to the next 4 items. If not, I couldn't explain why it does work :) – StepTNT May 14 '15 at 11:09
• @StepTNT Does it really work though? It didn't for me. Did you test it with test values against a working C version? – JS1 May 14 '15 at 11:24
• @StepTNT It made the explanation easier. I also changed the next line, which used to read movaps xmm7, [mask1 + edx*4 - 16] to movaps xmm7, [mask1 + edx - 16] (notice I don't multiply edx by 4 because I shifted it by two extra already). Both ways should be equivalent. – JS1 May 14 '15 at 17:21

This comment:

; ------------------------------------------------------------
; Entering the function
; ------------------------------------------------------------


is a little useless; it's it quite obvious that the start of the function/subroutine is right there.

However, the rest of your comments are very detailed and very well placed.

Don't use dec; using sub REG, 1 is faster.

dec runs in 1 CPU cycle. sub REG, 1 runs in half a CPU cycle.

Source: this post

Other than that, I'm not too familiar with xmm registers and the instructions that go with them.

However, as I look into them, if any improvements pop into my head, I'll update my post.

• I have that first comment because this one's part of a bigger file with multiple functions so I use them as separators between functions. Thanks for the link, I'll change all the inc/dec with add/sub. – StepTNT May 13 '15 at 21:12

If you are willing to address the masks through the movups instruction then a 50% gain in storage for these masks can be achieved. The reduced alignment condition will further add to this gain.

align 4
masks:  dd      0xffffffff, 0xffffffff, 0xffffffff, 0x00, 0x00, 0x00


Comparing EDX to 0 can also be done through decrementing EDX and watching for the sign flag. A simple negation of EDX gives the negative offset needed to fetch the correct mask.

.loopC:  dec    edx                     ;{0,1,2,3} -> {-1,0,1,2}

To clear the EBX register use xor ebx,ebx. It is shorter than mov ebx,0.
To jump if the ECX register is clear use jecxz .loopC. It is shorter than the combo cmp ecx,0 je .dU1. Please note that I changed the destination to cope with the possibility of the argument Dimensions being zero! Your code would erroneously have used mask3 in this case.
• Thak you for your comment, I like your idea on how to use masks! BTW there's no further check on input's size because it's guaranteed that there's at least one element, so Dimensions can't be 0 (I check for it in the caller function) – StepTNT May 19 '15 at 8:28