Fast median algorithm

[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

[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

Code Select Expand


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

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

Code Select Expand


/* 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];
}




[nidud]

deleted

[jj2007]

Here's one, (q)sort and grab the middle

Trivial but...

without 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

Code Select Expand


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 Expand


//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);



[daydreamer]

what about use MAXPS and MINPS instructions?

[nidud]

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.

[nidud]

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 ?

Code Select Expand



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];
}


[HSE]

Hi Guga!

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

[Siekmanski]

Hi guga,

Cooked up this algorithm.

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

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