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 "tpcextensions.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 doubleArrayMatch(
  const double *a1,
  const double *a2,
  const unsigned int n,
  const double lim
) {
  if(a1==NULL || a2==NULL || n<1) return 0;
  for(unsigned int i=0; i<n; i++) if(!doubleMatch(a1[i], a2[i], lim)) return 0;
  return 1;
}
/*****************************************************************************/
 
/*****************************************************************************/
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 doubleCopyFinite(
  double *t,
  double *s,
  const unsigned int n
) {
  unsigned int i, tn=0;
  if(t==NULL || s==NULL || n<1) return(0);
  for(i=0; i<n; i++) if(isfinite(s[i])) t[tn++]=s[i];
  return(tn);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int doubleNaNs(
  double *a,
  const unsigned int n
) {
  if(a==NULL || n<1) return(0);
  unsigned int i, nn=0;
  for(i=0; i<n; i++) if(isnan(a[i])) nn++;
  return(nn);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int doubleRange(
  double *a,
  const unsigned int n,
  double *amin,
  double *amax
) {
  if(amin!=NULL) *amin=nan("");
  if(amax!=NULL) *amax=nan("");
  if(a==NULL || n<1) return(0);
 
  double mi=nan(""), ma=nan("");
  unsigned int i, nn=0;
  for(i=0; i<n; i++) if(isfinite(a[i])) {
    nn++;
    if(!isfinite(mi) || mi>a[i]) mi=a[i];
    if(!isfinite(ma) || ma<a[i]) ma=a[i];
  }
  if(amin!=NULL) *amin=mi;
  if(amax!=NULL) *amax=ma;
  return(nn);
}
/*****************************************************************************/
 
/*****************************************************************************/
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);
}
/*****************************************************************************/
 
/*****************************************************************************/
double doubleWMean(
  double *a,
  double *w,
  const unsigned int n
) {
  if(a==NULL || w==NULL || n<1) return(nan(""));
  double s=0.0, ws=0.0;
  unsigned int i, ci=0;
  for(i=0; i<n; i++) if(isnormal(a[i]) && isnormal(w[i]) && w[i]>0.0) {
    ci++; s+=w[i]*a[i]; ws+=w[i];}
  if(ci<1 || !(ws>0.0)) return(nan(""));
  return(s/ws);
}
/*****************************************************************************/
 
/*****************************************************************************/
int doubleGetWithUnit(
  const char *s,
  double *v,
  int *u
) {
  if(v!=NULL) *v=nan(""); 
  if(u!=NULL) *u=UNIT_UNKNOWN; 
  if(s==NULL) return 1;
  int n=strTokenNr(s, " \t"); if(n==0 || n>2) return 1;
  if(n==1) return(atofCheck(s, v)); // no unit
  /* separate number and unit */
  char buf[64];
  n=strTokenNCpy(s, " \t", 1, buf, 64);
  if(atofCheck(buf, v)) return 1;
  if(u==NULL) return 0;
  n=strTokenNCpy(s, " \t", 2, buf, 64);
  *u=unitIdentify(buf);
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
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);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int doubleNonzeroes(
  double *a,
  const unsigned int n
) {
  if(a==NULL || n<1) return(0);
  unsigned int i, m=0;
  for(i=0; i<n; i++) if(fabs(a[i])>0.0) m++;
  return(m);
}
/*****************************************************************************/
 
/*****************************************************************************/
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);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int doubleAbsMaxIndex(
  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) || fabs(a[i])>mv) {mv=fabs(a[i]); mi=i;}
  return(mi);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int doubleMinIndex(
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int doubleAbsMinIndex(
  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) || fabs(a[i])<mv) {mv=fabs(a[i]); mi=i;}
  return(mi);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int doubleGEIndex(
  double *a,
  const unsigned int n,
  double lim
) {
  if(a==NULL || n<1) return(n);
  for(unsigned int i=0; i<n; i++) if(a[i]>=lim) return(i);
  return(n);
}
/*****************************************************************************/
unsigned int doubleGTIndex(
  double *a,
  const unsigned int n,
  double lim
) {
  if(a==NULL || n<1) return(n);
  for(unsigned int i=0; i<n; i++) if(a[i]>lim) return(i);
  return(n);
}
/*****************************************************************************/
 
/*****************************************************************************/
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));
}
/*****************************************************************************/
 
/*****************************************************************************/