Author Topic: Fast median algorithm  (Read 3318 times)

guga

  • Member
  • *****
  • Posts: 1274
  • Assembly is a state of art.
    • RosAsm
Fast median algorithm
« on: July 18, 2020, 02:22:35 AM »
Hi Guys

Someone knows a fast median algorithm for Integers and also for floats for sorted and unsorted values (without using sorting algorithms, btw) ? The inputed numbers are random and can be repeated. And the algorithm needs to adapt for even and odd lengths of any size.

Ex: A sequence of integers (and also floats/doubles): 1, 18, 0, 0, 0, 0, 2, 2, 5, 99, 600, 785, 758 , 17, 11, 21, 23, 24, 66, 21, 14, 57, 17

The goal is to obey this rules: https://www.mathsisfun.com/median.html

Coding in Assembly requires a mix of:
80% of brain, passion, intuition, creativity
10% of programming skills
10% of alcoholic levels in your blood.

My Code Sites:
http://rosasm.freeforums.org
http://winasm.tripod.com

jj2007

  • Member
  • *****
  • Posts: 10539
  • Assembler is fun ;-)
    • MasmBasic
Re: Fast median algorithm
« Reply #1 on: July 18, 2020, 02:54:31 AM »
You googled median algorithm and found Tony Hoare's algo, I suppose?

guga

  • Member
  • *****
  • Posts: 1274
  • Assembly is a state of art.
    • RosAsm
Re: Fast median algorithm
« Reply #2 on: July 18, 2020, 03:15:50 AM »
Hi JJ.

No i found this one.

https://stackoverflow.com/questions/810657/fastest-code-c-c-to-select-the-median-in-a-set-of-27-floating-point-values

Code: [Select]
// return the median value in a vector of 27 floats pointed to by a
float heapMedian3( float *a )
{
   float left[14], right[14], median, *p;
   unsigned char nLeft, nRight;

   // pick first value as median candidate
   p = a;
   median = *p++;
   nLeft = nRight = 1;

   for(;;)
   {
       // get next value
       float val = *p++;

       // if value is smaller than median, append to left heap
       if( val < median )
       {
           // move biggest value to the heap top
           unsigned char child = nLeft++, parent = (child - 1) / 2;
           while( parent && val > left[parent] )
           {
               left[child] = left[parent];
               child = parent;
               parent = (parent - 1) / 2;
           }
           left[child] = val;

           // if left heap is full
           if( nLeft == 14 )
           {
               // for each remaining value
               for( unsigned char nVal = 27 - (p - a); nVal; --nVal )
               {
                   // get next value
                   val = *p++;

                   // if value is to be inserted in the left heap
                   if( val < median )
                   {
                       child = left[2] > left[1] ? 2 : 1;
                       if( val >= left[child] )
                           median = val;
                       else
                       {
                           median = left[child];
                           parent = child;
                           child = parent*2 + 1;
                           while( child < 14 )
                           {
                               if( child < 13 && left[child+1] > left[child] )
                                   ++child;
                               if( val >= left[child] )
                                   break;
                               left[parent] = left[child];
                               parent = child;
                               child = parent*2 + 1;
                           }
                           left[parent] = val;
                       }
                   }
               }
               return median;
           }
       }

       // else append to right heap
       else
       {
           // move smallest value to the heap top
           unsigned char child = nRight++, parent = (child - 1) / 2;
           while( parent && val < right[parent] )
           {
               right[child] = right[parent];
               child = parent;
               parent = (parent - 1) / 2;
           }
           right[child] = val;

           // if right heap is full
           if( nRight == 14 )
           {
               // for each remaining value
               for( unsigned char nVal = 27 - (p - a); nVal; --nVal )
               {
                   // get next value
                   val = *p++;

                   // if value is to be inserted in the right heap
                   if( val > median )
                   {
                       child = right[2] < right[1] ? 2 : 1;
                       if( val <= right[child] )
                           median = val;
                       else
                       {
                           median = right[child];
                           parent = child;
                           child = parent*2 + 1;
                           while( child < 14 )
                           {
                               if( child < 13 && right[child+1] < right[child] )
                                   ++child;
                               if( val <= right[child] )
                                   break;
                               right[parent] = right[child];
                               parent = child;
                               child = parent*2 + 1;
                           }
                           right[parent] = val;
                       }
                   }
               }
               return median;
           }
       }
   }
}

But i didn´t tested it and also it seems to work only for 27 numbers, and also it seems to be kind of huge to be fast as intended to.
Coding in Assembly requires a mix of:
80% of brain, passion, intuition, creativity
10% of programming skills
10% of alcoholic levels in your blood.

My Code Sites:
http://rosasm.freeforums.org
http://winasm.tripod.com

guga

  • Member
  • *****
  • Posts: 1274
  • Assembly is a state of art.
    • RosAsm
Re: Fast median algorithm
« Reply #3 on: July 18, 2020, 03:54:35 AM »
I found also this oine that seems to be faster, but i´m unable to compile

http://blog.beamng.com/a-faster-selection-algorithm/
Code: [Select]
/* This source code is in the public domain */
#define F_SWAP(a,b) { float temp=(a);(a)=(b);(b)=temp; }
# Note: The code needs more than 2 elements to work
float lefselect(float a[], const int n, const int k) {
    int l=0, m = n-1, i=l, j=m;
    float x;
    while (l<m) {
        if( a[k] < a ) F_SWAP(a,a[k]);
        if( a[j] < a ) F_SWAP(a,a[j]);
        if( a[j] < a[k] ) F_SWAP(a[k],a[j]);
        x=a[k];
        while (j>k & i<k) {
            do i++; while (a<x);
            do j--; while (a[j]>x);
            F_SWAP(a,a[j]);
        }
        i++; j--;
        if (j<k) {
            while (a<x) i++;
            l=i; j=m;
        }
        if (k<i) {
            while (x<a[j]) j--;
            m=j; i=l;
        }
    }
    return a[k];
}

 
Coding in Assembly requires a mix of:
80% of brain, passion, intuition, creativity
10% of programming skills
10% of alcoholic levels in your blood.

My Code Sites:
http://rosasm.freeforums.org
http://winasm.tripod.com

nidud

  • Member
  • *****
  • Posts: 1980
    • https://github.com/nidud/asmc
Re: Fast median algorithm
« Reply #4 on: July 18, 2020, 04:40:16 AM »
Here's one, (q)sort and grab the middle.

Code: [Select]
;
; build: asmc64 -q -pe median.asm
;
include stdio.inc
include stdlib.inc
include tchar.inc

type typedef long_t ; item size

compare proto :ptr, :ptr {
    mov eax,[rcx]
    sub eax,[rdx]
    sbb ecx,ecx
    sbb ecx,-1
    and eax,eax
    cmovnz eax,ecx
    }

memxchg proto :ptr, :ptr {
    mov eax,[rcx]
    mov r8d,[rdx]
    mov [rdx],eax
    mov [rcx],r8d
    }

    .data
    array dd 1, 18, 0, 0, 0, 0, 2, 2, 5, 99, 600, 785, 758 , 17, 11, 21, 23, 24, 66, 21, 14, 57, 17

    .code

median proc uses rsi rdi rbx p:ptr, n:size_t

  local result:ptr

    mov rax,rdx

    .if rax > 1
       
        mov rsi,rcx
        dec rax
        mov ecx,type
        mul rcx
        lea rdi,[rsi+rax]
        xor r9d,r9d ; stack level

        lea rax,[rdi+rcx] ; middle from (hi - lo) / 2
        sub rax,rsi
        xor edx,edx
        div rcx
        shr rax,1
        mul rcx
        add rax,rsi
        mov result,rax

        .while 1
           
            mov rcx,type
            lea rax,[rdi+rcx]
            sub rax,rsi
            .ifnz
                xor edx,edx
                div rcx
                shr rax,1
                mul rcx
            .endif
           
            lea rbx,[rsi+rax]

            .ifsd compare(rsi, rbx) > 0
                memxchg(rsi, rbx)
            .endif
            .ifsd compare(rsi, rdi) > 0
                memxchg(rsi, rdi)
            .endif
            .ifsd compare(rbx, rdi) > 0
                memxchg(rbx, rdi)
            .endif

            mov r10,rsi
            mov r11,rdi

            .while 1

                add r10,type
                .if r10 < rdi

                    .continue .ifsd compare(r10, rbx) <= 0
                .endif

                .while 1

                    sub r11,type

                    .break .if r11 <= rbx
                    .break .ifsd compare(r11, rbx) <= 0
                .endw

                mov rcx,r11
                mov rax,r10
                .break .if rcx < rax
                memxchg(rcx, rax)

                .if rbx == r11

                    mov rbx,r10
                .endif
            .endw

            add r11,type

            .while 1

                sub r11,type

                .break .if r11 <= rsi
                .break .ifd compare(r11, rbx)
            .endw

            mov rdx,r10
            mov rax,r11
            sub rax,rsi
            mov rcx,rdi
            sub rcx,rdx

            .ifs rax < rcx

                mov rcx,r11

                .if rdx < rdi

                    push rdx
                    push rdi
                    inc r9d
                .endif

                .if rsi < rcx

                    mov rdi,rcx
                    .continue
                .endif
            .else
                mov rcx,r11

                .if rsi < rcx

                    push rsi
                    push rcx
                    inc r9d
                .endif

                .if rdx < rdi

                    mov rsi,rdx
                    .continue
                .endif
            .endif

            .break .if !r9d

            dec r9d
            pop rdi
            pop rsi
        .endw
        mov rcx,result
        mov eax,[rcx]
    .endif
    ret

median endp

main proc

    printf( "median: %d\n", median( &array, lengthof(array) ) )
    ret

main endp

    end _tstart

jj2007

  • Member
  • *****
  • Posts: 10539
  • Assembler is fun ;-)
    • MasmBasic
Re: Fast median algorithm
« Reply #5 on: July 18, 2020, 05:36:52 AM »
Here's one, (q)sort and grab the middle

Trivial but...
without using sorting algorithms

guga

  • Member
  • *****
  • Posts: 1274
  • Assembly is a state of art.
    • RosAsm
Re: Fast median algorithm
« Reply #6 on: July 18, 2020, 05:43:38 AM »
Tks Nidud

But, it is in 64 bits and also uses qsort. I´m looking for some that does not need to use sorting algorithms inside.

I found another one made by andy dansby that seems promissor and claims to be even faster then lefselect i posted previously, but i didn´t tested yet

https://www.mmbforums.com//viewtopic.php?f=22&t=4411&start=160

Code: [Select]
float Dansby_kth(float array[], const unsigned int length, const unsigned int kTHvalue)
{
#define F_SWAP(a,b) { float temp=(a);(a)=(b);(b)=temp; }

unsigned int left = 0;
unsigned int left2 = 0;
unsigned int right = length - 1;
unsigned int right2 = length - 1;
unsigned int kthPlus =  kTHvalue;
unsigned int kthMinus = kTHvalue;

while (right > left)
{
//Median of 3 optimization
if( array[right] < array[left] ) F_SWAP(array[left],array[right]);
if( array[right] < array[kTHvalue] ) F_SWAP(array[kTHvalue],array[right]);
if( array[kTHvalue] < array[left] ) F_SWAP(array[left],array[kTHvalue]);
//Median of 3 optimization

const float temp = array[kTHvalue];//this is pivot

kthPlus = kTHvalue + 1;
kthMinus = kTHvalue - 1;

while ((right2 > kTHvalue) && (left2 < kTHvalue))
{
do
{
left2++;
}while (array[left2] < temp);

do
{
right2--;
}while (array[right2] > temp);

F_SWAP(array[left2],array[right2]);
}
left2++;
right2--;
 
if (right2 < kthMinus)
{
while (array[left2] < temp)
{
left2++;
}
left = left2;
right2 = right;
kthMinus --;
}
else
if (kthPlus < left2)
{
while (array[right2] > temp)
{
right2--;
}
right = right2;
left2 = left;
kthPlus ++;
}
else
if( array[left] < array[right] )
{
F_SWAP(array[right],array[left]);
}



}

#undef F_SWAP
return array[kTHvalue];
}



call from

Code: [Select]
//points to the center value (median)
int middleElement = elementsinarray / 2;

Dansby_kth(redArray, elementsinarray,middleElement);
Dansby_kth(greenArray, elementsinarray,middleElement);
Dansby_kth(blueArray, elementsinarray,middleElement);

Coding in Assembly requires a mix of:
80% of brain, passion, intuition, creativity
10% of programming skills
10% of alcoholic levels in your blood.

My Code Sites:
http://rosasm.freeforums.org
http://winasm.tripod.com

daydreamer

  • Member
  • *****
  • Posts: 1350
  • building nextdoor
Re: Fast median algorithm
« Reply #7 on: July 18, 2020, 05:53:22 AM »
what about use MAXPS and MINPS instructions?
Quote from Flashdance
Nick  :  When you give up your dream, you die
*wears a flameproof asbestos suit*
Gone serverside programming p:  :D
I love assembly,because its legal to write
princess:lea eax,luke
:)

nidud

  • Member
  • *****
  • Posts: 1980
    • https://github.com/nidud/asmc
Re: Fast median algorithm
« Reply #8 on: July 18, 2020, 06:02:36 AM »
 :biggrin:

They are all based on sort algorithms (what else) so there's no way around that. You compare, swap and return the middle.

Quote
But, it is in 64 bits and also uses qsort. I´m looking for some that does not need to use sorting algorithms inside.

No, the algo do not use qsort.

guga

  • Member
  • *****
  • Posts: 1274
  • Assembly is a state of art.
    • RosAsm
Re: Fast median algorithm
« Reply #9 on: July 18, 2020, 09:02:57 AM »
Thank you Nidud, but do you have it for 32 Bits so i can test it ?

Hi Daydreamer. Using maxps/minps ? Sounds interesting. Not sure how to implement it yet, but i´ll give a try.
Coding in Assembly requires a mix of:
80% of brain, passion, intuition, creativity
10% of programming skills
10% of alcoholic levels in your blood.

My Code Sites:
http://rosasm.freeforums.org
http://winasm.tripod.com

nidud

  • Member
  • *****
  • Posts: 1980
    • https://github.com/nidud/asmc
Re: Fast median algorithm
« Reply #10 on: July 18, 2020, 10:17:58 PM »
From what I see the heapMedian3 sample have the advantage of not trashing the array but needs a local buffer. The rest of the samples are just variations of the same.

but do you have it for 32 Bits so i can test it ?

You may compile or past the C source to https://gcc.godbolt.org to get the output you need.

guga

  • Member
  • *****
  • Posts: 1274
  • Assembly is a state of art.
    • RosAsm
Re: Fast median algorithm
« Reply #11 on: July 19, 2020, 05:34:13 AM »
Hi Nidud


Ok, i´ll give a try and test the Heapmedian3 to see if it works as expected. However, as you said, it needs an extra buffer in memory which could be overkilling when we need to calculate the median of a larger array then only 27 bytes longs (say a array of floats, ints etc) with 100 Mb for example.

I[´ll try porting it to see if it works and also check for speed.

I was trying to port Dansby_kth that is very fast and don´t uses extra buffer, but it is producing a wrong result and i can´t fix that :(

Do you have any idea of what could be possibly wrong ?

Code: [Select]

float Dansby_kth(float array[], const unsigned int length, const unsigned int kTHvalue)
{
#define F_SWAP(a,b) { float temp=(a);(a)=(b);(b)=temp; }

unsigned int left = 0;
unsigned int left2 = 0;
unsigned int right = length - 1;
unsigned int right2 = length - 1;
unsigned int kthPlus = kTHvalue;
unsigned int kthMinus = kTHvalue;

while (right > left)
{
//Median of 3 optimization
if (array[right] < array[left]) F_SWAP(array[left], array[right]);
if (array[right] < array[kTHvalue]) F_SWAP(array[kTHvalue], array[right]);
if (array[kTHvalue] < array[left]) F_SWAP(array[left], array[kTHvalue]);
//Median of 3 optimization

const float temp = array[kTHvalue];//this is pivot

kthPlus = kTHvalue + 1;
kthMinus = kTHvalue - 1;

while ((right2 > kTHvalue) && (left2 < kTHvalue))
{
do
{
left2++;
} while (array[left2] < temp);

do
{
right2--;
} while (array[right2] > temp);

F_SWAP(array[left2], array[right2]);
}
left2++;
right2--;

if (right2 < kthMinus)
{
while (array[left2] < temp)
{
left2++;
}
left = left2;
right2 = right;
kthMinus--;
}
else
if (kthPlus < left2)
{
while (array[right2] > temp)
{
right2--;
}
right = right2;
left2 = left;
kthPlus++;
}
else
if (array[left] < array[right])
{
F_SWAP(array[right], array[left]);
}



}

#undef F_SWAP
return array[kTHvalue];
}
Coding in Assembly requires a mix of:
80% of brain, passion, intuition, creativity
10% of programming skills
10% of alcoholic levels in your blood.

My Code Sites:
http://rosasm.freeforums.org
http://winasm.tripod.com

HSE

  • Member
  • *****
  • Posts: 1377
  • <AMD>< 7-32>
Re: Fast median algorithm
« Reply #12 on: July 19, 2020, 06:08:51 AM »
Hi Guga!

Apparently the algorithm is for uneven data. Obviously not a coincidence, example have 27 observations.

Siekmanski

  • Member
  • *****
  • Posts: 2314
Re: Fast median algorithm
« Reply #13 on: July 19, 2020, 06:19:19 AM »
Hi guga,

Cooked up this algorithm.

New one: see http://masm32.com/board/index.php?topic=8671.msg94735#msg94735

« Last Edit: July 20, 2020, 08:55:13 AM by Siekmanski »
Creative coders use backward thinking techniques as a strategy.

jj2007

  • Member
  • *****
  • Posts: 10539
  • Assembler is fun ;-)
    • MasmBasic
Re: Fast median algorithm
« Reply #14 on: July 19, 2020, 06:50:13 AM »
Hi Guga,

Just for curiosity: what's the purpose of your question? Do you have a real use for this, or is it just for fun? For most applications, I guess sorting and taking the middle element would be fast enough...