Distance between two n-dimensional points (NASM)

Question

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     

        ; ------------------------------------------------------------
        ; Reading parameters from stack
        ; ------------------------------------------------------------
        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
        ; masked execution
        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
        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
        ; if not we jump to the square root
.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
        movaps  xmm7, [mask1]           
        jmp     .dU                     ; now that I have the correct mask I can jump to the actual execution

.dU2:   cmp     edx, 2
        jne     .dU3
        movaps  xmm7, [mask2]
        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 
        andps   xmm1, xmm7              ; masking
        andps   xmm2, xmm7              ; masking

        ; 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      
        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?

Answer 1

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

Your main loop:

.loopP: movaps  xmm1, [esi+ebx*4]
        movaps  xmm2, [edi+ebx*4]
        subps   xmm1, xmm2
        mulps   xmm1, xmm1
        dec     ecx
        haddps  xmm1, xmm1
        haddps  xmm1, xmm1
        addps   xmm0, xmm1
        jz      .loopC
        add     ebx, 4
        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
        add     ebx, 4
        mulps   xmm1, xmm1
        dec     ecx
        haddps  xmm1, xmm1
        haddps  xmm1, xmm1
        addps   xmm0, xmm1
        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.

Loading the mask

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

.dU1:   cmp     edx, 1          
        jne     .dU2
        movaps  xmm7, [mask1]           
        jmp     .dU

.dU2:   cmp     edx, 2
        jne     .dU3
        movaps  xmm7, [mask2]
        jmp     .dU

.dU3:   movaps  xmm7, [mask3]

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 */ };
        mask128 = *(__m128 *) &mask1[(remainder<<4)-16];

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.

Comments

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
        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.

Answer 2

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.

Answer 3

Defining the masks.

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

Now mask1 can be found at address [masks+8 -0]
Now mask2 can be found at address [masks+8 -4]
Now mask3 can be found at address [masks+8 -8].

Using the masks.

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}
         js     .endD                   ;No masked iterations
         neg    edx                     ;{0,1,2} -> {0,-1,-2}
         movups xmm7, [masks+8 + edx*4]

Size optimizations.

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.