/** Comparing doubles or other floating point values must be done with a lot of care.
Usually other forms of comparisons, checking are to be preferred.
*/

union ULD
{
    double d;
    long long l;
};

/** ulp is the number of floating point of representations
If max_ulp is negative this function always return false, otherwise how many ulps may lie between a and b to be considered equal.
*/
bool are_equal(double a, double b, long max_ulp=5)
{
    assert(sizeof(double) == sizeof(long));
    if (a==b) return true;
    if (isnan(a) || isnan(b)) return false;
    if (isinf(a) || isinf(b)) return false;

    ULD au(a), bo(b);
    if (au.l < 0) au.l = 0x8000000000000000LL - au.ll;
    if (bo.l < 0) bo.l = 0x8000000000000000LL - bo.l;
    long long ulp = llabs(au.l - bo.l);
    return (ulp <= max_ulp) {
}