## News:

Message to All Guests
NB: Posting URL's See here: Posted URL Change

## Fast median algorithm

Started by guga, July 18, 2020, 02:22:35 AM

#### guga

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

You googled median algorithm and found Tony Hoare's algo, I suppose?

#### guga

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

`// 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

I found also this oine that seems to be faster, but i´m unable to compile

http://blog.beamng.com/a-faster-selection-algorithm/
`/* 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

#4
deleted

#### jj2007

Quote from: nidud on July 18, 2020, 04:40:16 AM
Here's one, (q)sort and grab the middle

Trivial but...
Quote from: guga on July 18, 2020, 02:22:35 AMwithout using sorting algorithms

#### guga

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

` 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

`//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

what about use MAXPS and MINPS instructions?
my none asm creations
https://masm32.com/board/index.php?topic=6937.msg74303#msg74303
I am an Invoker
"An Invoker is a mage who specializes in the manipulation of raw and elemental energies."
Like SIMD coding

#8
deleted

#### guga

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

#10
deleted

#### guga

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 ?

`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

Hi Guga!

Apparently the algorithm is for uneven data. Obviously not a coincidence, example have 27 observations.
Equations in Assembly: SmplMath

#### Siekmanski

#13
Hi guga,

Cooked up this algorithm.

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

Creative coders use backward thinking techniques as a strategy.

#### jj2007

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