// equalish.h : Functions for comparing floating-point numbers generically.
//

#ifndef __EQUALISH_H__
#define __EQUALISH_H__

// Inspired by the discussion found at:
// http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

/////////////////////////////////////////////////////////////////////////////

#define equalish AlmostEqual2sComplement

/////////////////////////////////////////////////////////////////////////////

bool AlmostEqualUlps(float A, float B, int maxUlps=10);
bool AlmostEqual2sComplement(float A, float B, int maxUlps=10);

#ifdef _DEBUG

bool test_equalish();

#endif // _DEBUG

/////////////////////////////////////////////////////////////////////////////

#endif // __EQUALISH_H__

//------------------------------------------------>8--( cut here )-----------
// equalish.cpp

#include "equalish.h"

bool AlmostEqualUlps(float A, float B, int maxUlps /* = 10 */)
{
    assert(sizeof(float) == sizeof(int));
    if (A == B)
        return true;
    int intDiff = abs(*(int*)&A - *(int*)&B);
    if (intDiff <= maxUlps)
        return true;
    return false;
}

bool AlmostEqual2sComplement(float A, float B, int maxUlps /* = 10 */)
{
    // Make sure maxUlps is non-negative and small enough that the
    // default NAN won't compare as equal to anything.
    assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);
    int aInt = *(int*)&A;
    // Make aInt lexicographically ordered as a twos-complement int
    if (aInt < 0)
        aInt = 0x80000000 - aInt;
    // Make bInt lexicographically ordered as a twos-complement int
    int bInt = *(int*)&B;
    if (bInt < 0)
        bInt = 0x80000000 - bInt;
    int intDiff = abs(aInt - bInt);
    if (intDiff <= maxUlps)
        return true;
    return false;
}

/////////////////////////////////////////////////////////////////////////////

#ifdef _DEBUG

#include <stdio.h>

bool test_equalish()
{
        bool success = TRUE;

        try
        {
                int i;
                float f;
                for (i = 0, f = -1.0; f <= 1.0; f += 1.0/9.0, i++)
                {
                        // f is far enough from zero that tiny offsets are insignificant
                        assert(equalish(f - 0.000000000000001, f + 0.000000000000001));
                        if (i == 9) 
                        {
                                // f should be close to zero; tiny offsets are more significant
                                assert(!equalish(f - 0.000000000001, f + 0.0000000000001));
                        }
                }
        }
        catch(exception e)
        {
                printf("exception: %s\n", e.what());
                success = FALSE;
        }

        return success;
}

#endif // _DEBUG