There is a build in MATLAB implementation for float difference in ULP ?
u = abs(A-B) / (abs(B) * eps(B))
u = (y - Y)/eps(abs(Y))
other ulps for transcendental by Intel
error = abs( (f(x)- ref_f(x)) / 2^(k-63) ) where 1<= 2^(-k)*f(x) < 2
Intel 64 and IA-32 Architectures Software
8.3.10 Transcendental Instruction Accuracy
boost C++ lib
Function float_distance finds the number of gaps/bits/ULP between any two floating-point values.
If the significands of floating-point numbers are viewed as integers, then their difference is the number of ULP/gaps/bits different.
Description - float_distance
Returns the distance between a and b: the result is always a signed integer value (stored in floating-point type FPT) representing the number of distinct representations between a and b.
Note that
float_distance(a, a) always returns 0.
float_distance(float_next(a), a) always returns 1.
float_distance(float_prior(a), a) always returns -1.
The function float_distance is equivalent to calculating the number of ULP (Units in the Last Place) between a and b except that it returns a signed value indicating whether a > b or not.
If the distance is too great then it may not be able to be represented as an exact integer by type FPT, but in practice this is unlikely to be a issue.
Another implementation based on find direction than search identity by counting ULP steps:
template <typename T>
GLM_FUNC_QUALIFIER uint float_distance(T const & x, T const & y)
{
uint ulp = 0;
if(x < y)
{
T temp = x;
while(temp != y)// && ulp < std::numeric_limits<std::size_t>::max())
{
++ulp;
temp = next_float(temp);
}
}
else if(y < x)
{
T temp = y;
while(temp != x)// && ulp < std::numeric_limits<std::size_t>::max())
{
++ulp;
temp = next_float(temp);
}
}
else // ==
{
}
return ulp;
}
bool relativelyEqual(float a, float b,
float maxRelativeDiff = FLT_EPSILON)
{
const float difference = fabs(a - b);
// Scale to the largest value.
a = fabs(a);
b = fabs(b);
const float scaledEpsilon =
maxRelativeDiff * max(a, b);
return difference <= scaledEpsilon;
}
float relative_difference(float a, float b)
{
return fabs((a - b) / min(a, b));
}
float epsilon_difference(float a, float b)
{
return relative_difference(a, b) /
FLT_EPSILON;
}
int32_t ulpsDistance(const float a, const float b)
{
// Save work if the floats are equal.
// Also handles +0 == -0
if (a == b) return 0;
const auto max =
std::numeric_limits<int32_t>::max();
// Max distance for NaN
if (isnan(a) || isnan(b)) return max;
// If one's infinite and they're not equal, max distance.
if (isinf(a) || isinf(b)) return max;
int32_t ia, ib;
memcpy(&ia, &a, sizeof(float));
memcpy(&ib, &b, sizeof(float));
// Don't compare differently-signed floats.
if ((ia < 0) != (ib < 0)) return max;
// Return the absolute value of the distance in ULPs.
int32_t distance = ia - ib;
if (distance < 0) distance = -distance;
return distance;
}
public static float ulp(float f)
Returns the size of an ulp of the argument.
An ulp, unit in the last place, of a float value is the positive distance between this floating-point value and the float value next larger in magnitude.
Note that for non-NaN x, ulp(-x) == ulp(x).
Special Cases:
If the argument is NaN, then the result is NaN.
If the argument is positive or negative infinity, then the result is positive infinity.
If the argument is positive or negative zero, then the result is Float.MIN_VALUE.
If the argument is ±Float.MAX_VALUE, then the result is equal to 2104.
Parameters:
f - the floating-point value whose ulp is to be returned
Returns:
the size of an ulp of the argument
Since:
1.5
bool AlmostEqualUlpsAndAbs(float A, float B,
float maxDiff, int maxUlpsDiff)
{
// Check if the numbers are really close -- needed
// when comparing numbers near zero.
float absDiff = fabs(A - B);
if (absDiff <= maxDiff)
return true;
Float_t uA(A);
Float_t uB(B);
// Different signs means they do not match.
if (uA.Negative() != uB.Negative())
return false;
// Find the difference in ULPs.
int ulpsDiff = abs(uA.i - uB.i);
if (ulpsDiff <= maxUlpsDiff)
return true;
return false;
}
Other references:
Knuth D.E. The art of computer programming, vol II, section 4.2, especially Floating-Point Comparison 4.2.2, pages 198-220.
David Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic"
Alberto Squassabia, Comparing floats listing
Google Floating-Point_Comparison guide
Boost.Test Floating-Point_Comparison
