doubleutil.c File Reference

Working with doubles. More...

#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"

Go to the source code of this file.

Functions
int	doubleMatch (const double v1, const double v2, const double lim)
int	doubleMatchRel (const double v1, const double v2, const double lim)
double	doubleMachEps ()
void	doubleCopy (double *t, double *s, const unsigned int n)
unsigned int	doubleMaxIndex (double *a, const unsigned int n)
double	doubleSum (double *a, const unsigned int n)
double	doubleMean (double *a, const unsigned int n)
int	doubleSpanPositives (double *a, const int n)
int	doubleCSpanPositives (double *a, const int n)
double	inverfc (double x)
void	statSortDouble (double *data, unsigned int n, int order)
void	statSortFloat (float *data, unsigned int n, int order)

Detailed Description

Working with doubles.

Author

Vesa Oikonen

Definition in file doubleutil.c.

Function Documentation

doubleCopy()

void doubleCopy	(	double *	t,
		double *	s,
		const unsigned int	n )

Copy double values from the 2nd array to the first.

Parameters

t	Target array
s	Source array
n	Length of arrays

Definition at line 93 of file doubleutil.c.

  {
  unsigned int i;
  if(t==NULL || s==NULL || n<1) return;
  for(i=0; i<n; i++) t[i]=s[i];
}

doubleCSpanPositives()

int doubleCSpanPositives	(	double *	a,
		const int	n )

Returns the length of array consisting of non-positive (<=0 and NaN) values.

Returns

Returns the index of first positive value in the given array.

See also

doubleSpanPositives

Parameters

a	Pointer to the array.
n	Length of the array.

Definition at line 186 of file doubleutil.c.

  {
  if(a==NULL || n<1) return(0);
  int i=0;
  for(i=0; i<n; i++) if(a[i]>0.0) break;
  return(i);
}

doubleMachEps()

double doubleMachEps

(

)

Estimates the machine epsilon, the upper bound on the relative error due to rounding in floating point arithmetic, within one order of magnitude of the true machine epsilon.

Standard C library should also have DBL_EPSILON in float.h.

Returns

Estimate of machine epsilon.

Author

Vesa Oikonen

See also

doubleMatchRel, doubleMatch

Definition at line 83 of file doubleutil.c.

{
  double macheps=1.0;
  do {macheps/=2.0;} while((1.0+macheps/2.0)!=1.0);
  return(macheps);
}

Referenced by householder_transform().

doubleMatch()

int doubleMatch	(	const double	v1,
		const double	v2,
		const double	lim )

Verifies that given two doubles have the same value inside given limits.

Values are considered to match also if both are NaNs.

Returns

1 if match is found, and 0 if not.

Author

Vesa Oikonen

See also

doubleMatchRel, doubleMachEps

Parameters

v1	First value
v2	Second value
lim	Limit for absolute difference (if <0 then test will fail every time)

Definition at line 28 of file doubleutil.c.

  {
  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;
}

doubleMatchRel()

int doubleMatchRel	(	const double	v1,
		const double	v2,
		const double	lim )

Verifies that given two doubles have the same value inside given relative limit |2*(v1-v2)/(v1+v2)|.

Values are considered to match also if both are NaNs, but not if mean is zero (unless both are exactly zero).

Returns

1 if match is found, and 0 if not.

Author

Vesa Oikonen

See also

doubleMatch, doubleMachEps

Parameters

v1	First value
v2	Second value
lim	Limit for relative difference (if <0 then test will fail every time)

Definition at line 55 of file doubleutil.c.

  {
  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;
}

doubleMaxIndex()

unsigned int doubleMaxIndex	(	double *	a,
		const unsigned int	n )

Find the maximum value in given double array.

Returns

The index [0..n-1] of maximum value in array; 0 is returned also in case of errors.

Parameters

a	Pointer to double array.
n	Length of array.

Definition at line 111 of file doubleutil.c.

  {
  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);
}

doubleMean()

double doubleMean	(	double *	a,
		const unsigned int	n )

Calculate the mean of values in given double array.

Returns

The mean of array values, or NaN in case of an error.

See also

doubleSum

Parameters

a	Pointer to double array.
n	Length of array.

Definition at line 148 of file doubleutil.c.

  {
  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);
}

doubleSpanPositives()

int doubleSpanPositives	(	double *	a,
		const int	n )

Returns the length of array consisting of only positive (>0 and not NaN) values.

Returns

Returns the index of first nonpositive value in the given array.

See also

doubleCSpanPositives

Parameters

a	Pointer to the array.
n	Length of the array.

Definition at line 168 of file doubleutil.c.

  {
  if(a==NULL || n<1) return(0);
  int i=0;
  for(i=0; i<n; i++) if(!(a[i]>0.0)) break;
  return(i);
}

doubleSum()

double doubleSum	(	double *	a,
		const unsigned int	n )

Calculate the sum of values in given double array.

Returns

The sum of array values.

See also

doubleMean

Parameters

a	Pointer to double array.
n	Length of array.

Definition at line 130 of file doubleutil.c.

  {
  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);
}

inverfc()

double inverfc

(

double

x

)

Inverse complementary error function erfc^{-1}(x).

Based on the code by Acklam PJ, 2010.

Returns

Returns an approximate value of erfc^{-1}(x).

Parameters

x	Parameter for the inverse complementary error function.

Definition at line 205 of file doubleutil.c.

  {
  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));
}

Referenced by imgGaussianFIRFilter().

statSortDouble()

void statSortDouble	(	double *	data,
		unsigned int	n,
		int	order )

Sort the given double array into ascending or descending order.

Author

Vesa Oikonen

Parameters

data	Pointer to data array of size n
n	Length of data array
order	Ascending (0) or descending (<>0) order

Definition at line 267 of file doubleutil.c.

  {
  if(n<2 || data==NULL) return;
  if(order==0) qsort(data, n, sizeof(double), statDoubleCompAsc);
  else qsort(data, n, sizeof(float), statDoubleCompDesc);
}

statSortFloat()

void statSortFloat	(	float *	data,
		unsigned int	n,
		int	order )

Sort the given float array into ascending or descending order.

Author

Vesa Oikonen

Parameters

data	Pointer to data array of size n
n	Length of data array
order	Ascending (0) or descending (<>0) order

Definition at line 301 of file doubleutil.c.

  {
  if(n<2 || data==NULL) return;
  if(order==0) qsort(data, n, sizeof(float), statFloatCompAsc);
  else qsort(data, n, sizeof(float), statFloatCompDesc);
}