00001 #include "mrilib.h"
00002 
00003 /*------------------------------------------------------------------------
00004    Compute the median of an array of floats.  Will rearrange (partially
00005    sort) the array in the process.  The algorithm is based on Quicksort,
00006    where we only keep the partition that has the middle element.
00007    For large n, this is faster than finishing the whole sort.
00008 --------------------------------------------------------------------------*/
00009 
00010 /* macro for median-of-3 */
00011 
00012 #undef  MED3
00013 #define MED3(a,b,c) ( (a < b) ? ( (a > c) ? a : (b > c) ? c : b )  \
00014                               : ( (b > c) ? b : (a > c) ? c : a ) )
00015 
00016 /* macro for swap (duh) */
00017 
00018 #undef  SWAP
00019 #define SWAP(x,y) (temp=(x),(x)=(y),(y)=temp)
00020 
00021 float qmed_float( int n , float * ar )
00022 {
00023    register int i , j ;           /* scanning indices */
00024    register float temp , pivot ;  /* holding places */
00025    register float * a = ar ;
00026 
00027    int left , right , mid , nodd ;
00028 
00029    switch( n ){
00030       case 1: return ar[0] ;
00031       case 2: return 0.5*(ar[0]+ar[1]) ;
00032       case 3: return MED3( ar[0] , ar[1] , ar[2] ) ;
00033    }
00034 
00035    left = 0 ; right = n-1 ; mid = n/2 ; nodd = ((n & 1) != 0) ;
00036 
00037    /* loop while the subarray is at least 3 long */
00038 
00039    while( right-left > 1  ){  /* work on subarray from left -> right */
00040 
00041       i = ( left + right ) / 2 ;   /* middle of subarray */
00042 
00043       /* sort the left, middle, and right a[]'s */
00044 
00045       if( a[left] > a[i]     ) SWAP( a[left]  , a[i]     ) ;
00046       if( a[left] > a[right] ) SWAP( a[left]  , a[right] ) ;
00047       if( a[i] > a[right]    ) SWAP( a[right] , a[i]     ) ;
00048 
00049       pivot = a[i] ;                 /* a[i] is the median-of-3 pivot! */
00050       a[i]  = a[right] ;
00051 
00052       i = left ;                     /* initialize scanning */
00053       j = right ;
00054 
00055       /*----- partition:  move elements bigger than pivot up and elements
00056                           smaller than pivot down, scanning in from ends -----*/
00057 
00058       do{
00059         for( ; a[++i] < pivot ; ) ;  /* scan i up,   until a[i] >= pivot */
00060         for( ; a[--j] > pivot ; ) ;  /* scan j down, until a[j] <= pivot */
00061 
00062         if( j <= i ) break ;         /* if j meets i, quit */
00063 
00064         SWAP( a[i] , a[j] ) ;
00065       } while( 1 ) ;
00066 
00067       /*----- at this point, the array is partitioned -----*/
00068 
00069       a[right] = a[i] ;           /* restore the pivot */
00070       a[i]     = pivot ;
00071 
00072       if( i == mid ){             /* good luck */
00073          if( nodd ) return pivot ;   /* exact middle of array */
00074 
00075          temp = a[left] ;            /* must find next largest element below */
00076          for( j=left+1 ; j < i ; j++ )
00077             if( a[j] > temp ) temp = a[j] ;
00078 
00079          return 0.5*(pivot+temp) ;   /* return average */
00080       }
00081 
00082       if( i <  mid ) left  = i ; /* throw away bottom partition */
00083       else           right = i ; /* throw away top partition    */
00084 
00085    }  /* end of while sub-array is long */
00086 
00087    return (nodd) ? a[mid] : 0.5*(a[mid]+a[mid-1]) ;
00088 }
00089 
00090 int main( int argc , char * argv[] )
00091 {
00092    int nn,ii , med , mrep=1 , jj ;
00093    float * ar , * br ;
00094    float mm ;
00095    double ct , ctsum=0.0 ;
00096 
00097    if( argc < 2 ){printf("Usage: qmed N [m]\n");exit(0);}
00098    nn = strtol(argv[1],NULL,10) ; if( nn == 0 ) exit(1);
00099 
00100    if( argc == 3 ){
00101       mrep = strtol(argv[2],NULL,10) ; if( mrep < 1 ) mrep = 1 ;
00102    }
00103 
00104    med = ( nn > 0 ) ; if( !med ) nn = -nn ;
00105 
00106    ar = (float *) malloc( sizeof(float)*nn ) ;
00107    br = (float *) malloc( sizeof(float)*nn ) ;
00108    srand48((long)(237+jj)) ;
00109    for( ii=0 ; ii < nn ; ii++ ) br[ii] = (float)(lrand48() % nn)*0.1 ;
00110 
00111    ct = COX_cpu_time() ;
00112    for( jj=0 ; jj < mrep ; jj++ ){
00113 
00114       memcpy( ar , br , sizeof(float)*nn ) ;
00115 
00116       if( med ){
00117          mm = qmed_float( nn , ar ) ;
00118       } else {
00119          qsort_float( nn , ar ) ;
00120          if( nn%2 == 1 ) mm = ar[nn/2] ;
00121          else            mm = 0.5*(ar[nn/2]+ar[nn/2-1]) ;
00122       }
00123    }
00124    ct = COX_cpu_time() - ct ;
00125    printf("Median = %g  CPU = %g\n",mm,ct) ;
00126    exit(0) ;
00127 }