### Author Topic: Fast median algorithm  (Read 3318 times)

#### guga

• Member
• Posts: 1274
• Assembly is a state of art.
##### 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 ;-)
##### 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.
##### 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 afloat 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.
##### 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 workfloat 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
##### 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.incinclude stdlib.incinclude tchar.inctype typedef long_t ; item sizecompare 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    .codemedian 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    retmedian endpmain proc    printf( "median: %d\n", median( &array, lengthof(array) ) )    retmain endp    end _tstart`

#### jj2007

• Member
• Posts: 10539
• Assembler is fun ;-)
##### 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.
##### 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
##### Re: Fast median algorithm
« Reply #8 on: July 18, 2020, 06:02:36 AM »

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.
##### 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
##### 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.
##### 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 ;-)
##### 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...