doubleutil.c

Go to the documentation of this file.

 
/*****************************************************************************/
 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <ctype.h>
#include <unistd.h>
#include <math.h>
#include "libtpcmisc.h"
/*****************************************************************************/
 
/*****************************************************************************/
int doubleMatch(
  const double v1,
  const double v2,
  const double lim
) {
  if(isnan(v1) && isnan(v2)) return 1;
  if(isnan(v1) || isnan(v2)) return 0;
  if(v1==v2) return 1;
  if(isnan(lim) || lim<0.0) return 0;
  if(fabs(v1-v2)<=lim) return 1;
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int doubleMatchRel(
  const double v1,
  const double v2,
  const double lim
) {
  if(isnan(v1) && isnan(v2)) return 1;
  if(isnan(v1) || isnan(v2)) return 0;
  if(v1==v2) return 1;
  if(isnan(lim)) return 0;
  double mean;
  mean=0.5*(v1+v2); if(!isnormal(mean)) return 0;
  if(fabs((v1-v2)/mean)<=lim) return 1;
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
double doubleMachEps()
{
  double macheps=1.0;
  do {macheps/=2.0;} while((1.0+macheps/2.0)!=1.0);
  return(macheps);
}
/*****************************************************************************/
 
/*****************************************************************************/
void doubleCopy(
  double *t,
  double *s,
  const unsigned int n
) {
  unsigned int i;
  if(t==NULL || s==NULL || n<1) return;
  for(i=0; i<n; i++) t[i]=s[i];
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int doubleMaxIndex(
  double *a,
  const unsigned int n
) {
  if(a==NULL || n<1) return(0);
  unsigned int i, mi=0;
  double mv=nan("");
  for(i=0; i<n; i++) if(isnan(mv) || a[i]>mv) {mv=a[i]; mi=i;}
  return(mi);
}
/*****************************************************************************/
 
/*****************************************************************************/
double doubleSum(
  double *a,
  const unsigned int n
) {
  double s=0.0;
  if(a==NULL || n<1) return(s);
  for(unsigned int i=0; i<n; i++) if(!isnan(a[i])) s+=a[i];
  return(s);
}
/*****************************************************************************/
 
/*****************************************************************************/
double doubleMean(
  double *a,
  const unsigned int n
) {
  if(a==NULL || n<1) return(nan(""));
  double s=0.0;
  unsigned int i, ci=0;
  for(i=0; i<n; i++) if(!isnan(a[i])) {ci++; s+=a[i];}
  if(ci<1) return(nan(""));
  return(s/(double)ci);
}
/*****************************************************************************/
 
/*****************************************************************************/
int doubleSpanPositives(
  double *a,
  const int n
) {
  if(a==NULL || n<1) return(0);
  int i=0;
  for(i=0; i<n; i++) if(!(a[i]>0.0)) break;
  return(i);
}
/*****************************************************************************/
 
/*****************************************************************************/
int doubleCSpanPositives(
  double *a,
  const int n
) {
  if(a==NULL || n<1) return(0);
  int i=0;
  for(i=0; i<n; i++) if(a[i]>0.0) break;
  return(i);
}
/*****************************************************************************/
 
/*****************************************************************************/
double inverfc(
  double x
) {
  static const double a[] = {
    -3.969683028665376E+01, 2.209460984245205E+02, -2.759285104469687E+02,
    1.383577518672690E+02, -3.066479806614716E+01,  2.506628277459239
  };
  static const double b[] = {
    -5.447609879822406E+01, 1.615858368580409E+02, -1.556989798598866E+02, 
    6.680131188771972E+01, -1.328068155288572E+01
  };
  static const double c[] = {
    -7.784894002430293E-03, -3.223964580411365E-01, -2.400758277161838, 
    -2.549732539343734, 4.374664141464968, 2.938163982698783
  };
  static const double d[] = {
    7.784695709041462E-03, 3.224671290700398E-01, 2.445134137142996, 3.754408661907416
  };
  double y, e, u;
  x/=2.0;
  if(x<0.02425) {
    double q=sqrt(-2.0 * log(x));
    y = (((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5])
         / ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0);
  } else if(x<=0.97575) {
    double q=x-0.5;
    double r=q*q;
    y = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q
         / (((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0);
  } else {
    double q=sqrt(-2.0 * log(1.0-x));
    y = -(((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5])
          / ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0);
  }
  e=0.5*erfc(-y/M_SQRT2)-x;
  u=(e)*M_SQRT2*M_PI*exp(0.5*y*y);
  y-=u/(1.0 + 0.5*y*u);
  return(fabs(y/M_SQRT2));
}
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
static int statDoubleCompAsc(const void *i, const void *j)
{
  const double *di = (const double *)i;
  const double *dj = (const double *)j;
  return(*di > *dj) - (*di < *dj);
}
static int statDoubleCompDesc(const void *i, const void *j)
{
  const double *di = (const double *)i;
  const double *dj = (const double *)j;
  return(*di < *dj) - (*di > *dj);
}
 
void statSortDouble(
  double *data,
  unsigned int n,
  int order
) {
  if(n<2 || data==NULL) return;
  if(order==0) qsort(data, n, sizeof(double), statDoubleCompAsc);
  else qsort(data, n, sizeof(float), statDoubleCompDesc);
}
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
static int statFloatCompAsc(const void *i, const void *j)
{
  const float *di = (const float *)i;
  const float *dj = (const float *)j;
  return(*di > *dj) - (*di < *dj);
}
static int statFloatCompDesc(const void *i, const void *j)
{
  const float *di = (const float *)i;
  const float *dj = (const float *)j;
  return(*di < *dj) - (*di > *dj);
}
 
void statSortFloat(
  float *data,
  unsigned int n,
  int order
) {
  if(n<2 || data==NULL) return;
  if(order==0) qsort(data, n, sizeof(float), statFloatCompAsc);
  else qsort(data, n, sizeof(float), statFloatCompDesc);
}
/*****************************************************************************/
 
/*****************************************************************************/