News:

Masm32 SDK description, downloads and other helpful links
Message to All Guests

Main Menu

Fast median algorithm

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

Previous topic - Next topic

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

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

#4
deleted

jj2007


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

nidud

#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

nidud

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