#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_errno.h>


#include <math.h>


/*
 * Recursively calculates the mean of n doubles according to the
 * folowing equation:
 *   mean(n) = (mean(n-1) * (n-1) + a(n)) / n
 */
//double mean_double2(double* data, int startpos,int length ){
                 /*  double* data  data to arrange                */
                 /*  int startpos  one-based start index      	  */
                 /*  int length    length of the subarray to sort */
//  int i, k;
//  double mean;
//  for(i = startpos, mean = data[startpos - 1], k = 0; i < length; i++, k++) {
//    mean = (i * mean + data[k]) / (i + 1);
//  }  
//  return mean;
//}

 /* Calculates standart deviation */
/* double sd_double2( double *data, double data_mean, int n_elements ){
  int i;
  double *x;
  x=malloc(n_elements*sizeof(double));
  for(i=0;i<n_elements;i++){
   x[i]=(data[i]-data_mean)*(data[i]-data_mean);
  }
  return sqrt(mean_double(x,0,n_elements)); 
 }*/


/*
 * Comparer for QSort.
 */
#define COMPARER_EPS 0.000001 /* accuracy. Probably, shoul be tuned */
/*__compar_fn_t compare_double(double *a, double *b) {                               
 int res1=1,resm1=-1,res0=0;
 
 if (((*a)-(*b))>(COMPARER_EPS)){
  return &res1;
 }                    
 else {
  if( ((*a)-(*b))<-1*(COMPARER_EPS) ){
   return &resm1;
  }
  else{       
   return &res0;
  }
 }
} */                                                      
/*
 * Arranges doubles within the array
 */
void sort_data_double2(double* data, int startpos, int length) {
 //--tmp--
  int i;
/*  for(i=0;i<length;i++){
   fprintf(stderr,"%lf\n",data[i]);
  }*/
 //-------
 int j;
 double tmp;
 for(i=0;i<length;i++){
  for(j=0;j<length;j++){
   if( data[i]<data[j] ){
    tmp=data[j];
    data[j]=data[i];
    data[i]=tmp;
   }
  }
 }
 //qsort(data, sizeof(double), length, &compare_double);
 //---
 /*printf("--------------\n");
   for(i=0;i<length;i++){
      fprintf(stderr,"%lf\n",data[i]);
        }*/
         //-------
         
 //---
 
 
}
                                                              
/*
 * Calculates the median of the SORTED data array
 * (use sort_data_double())
 */
double median_double2(double* data,
                     int startpos,
                     int length
                     ) {

 sort_data_double2( data, startpos, length);

 if (length % 2) {
  /* FIXME: too inevident... */
//  return (data[startpos+length/2]+data[startpos+length/2-1])/2;
  return data[length/2];
 }
 else {
  /* FIXME: should work fine?? */
//  return data[startpos+length/2];
  return (data[(int)(length/2-0.5)]+data[(int)(length/2+0.5)])/2;
 }
}

double median_sigma_m2(double mean_X,double *X,int N){
 double Xm,res;
 double XXX[10000];
 int i;
 for(i=0;i<N;i++){
  XXX[i]=fabs(X[i]-mean_X);
 }
 Xm=gsl_stats_mean(XXX,1,N);
 res=Xm/(0.674*sqrt(1+M_PI/(2*N)));
// fprintf(stderr,"median_sigma=%lf median=%lf mean=%lf\n",res,Xm,mean_X);
 return res;
}