median.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
/*****************************************************************************/
#include "tpcstatist.h"
/*****************************************************************************/
 
/*****************************************************************************/
double statMedian(
  double *a,
  const int n
) {
  if(a==NULL || n<1) return(nan(""));
  if(n==1) return(a[0]);
 
  if(n%2) {
    int k=(n-1)/2; 
    return(statKthSmallest(a, n, k));
  } else {
    int k=n/2; 
    double d1=statKthSmallest(a, n, k-1);
    double d2=statKthSmallest(a, n, k);
    return(0.5*(d1+d2));
  }
}
/*****************************************************************************/
 
/*****************************************************************************/
double statKthSmallest(
  double *a,
  const int n,
  const int k
) {
  //printf("\n%s(a[], %u, %u)\n", __func__, n, k); fflush(stdout);
  if(a==NULL || n<1 || k<0 || k>=n) return(nan(""));
  if(n==1) return(a[0]);
 
  int l=0, m=n-1; // note that algorithm does not work with unsigned integers
  while(l<m) {
    int i=l, j=m;
    double x=a[k]; 
    do {
      while(a[i]<x) i++;
      while(x<a[j]) j--;
      if(i<=j) {double s=a[i]; a[i]=a[j]; a[j]=s; i++; j--;}
    } while(i<=j);
    if(j<k) l=i;
    if(k<i) m=j;
  }
  return(a[k]);
}
/*****************************************************************************/
 
/*****************************************************************************/