#include <stdlib.h>
#include <stdio.h>
#include <fenv.h>
// compile with gcc -std=c99 for fpclassify
inline int classify(double x) // filter out bad numbers correct for inexact, INT_MIN is the error return;
{
switch(fpclassify(x))
{
case FP_INFINITE:
errno=ERANGE;
return INT_MIN;
case FP_NAN:
errno=EINVAL;
return INT_MIN;
case FP_NORMAL:
return 2;
case FP_SUBNORMAL:
case FP_ZERO:
return 0;
}
return FP_NORMAL; //default
}
// doug gwyn's reldif function
// usage: if(reldif(a, b) <= TOLERANCE) ...
#define Abs(x) ((x) < 0 ? -(x) : (x))
#define Max(a, b) ((a) > (b) ? (a) : (b))
double maxulp=0; // tweek this into TOLERANCE
inline double reldif(double a, double b)
{
double c = Abs(a);
double d = Abs(b);
int rc=0;
d = Max(c, d);
d = (d == 0.0) ? 0.0 : Abs(a - b) / d;
rc=classify(d);
if(!rc) // correct almost zeroes to zero
d=0.;
if(rc == INT_MIN )
{ // error return
errno=ERANGE;
d=DBL_MAX;
perror("Error comparing values");
}
return d;
}
// usage: if(reldif(a, b) <= TOLERANCE) ...