Page 40 - Developer

INNER PRODUCT // bRUCE DawsON

listing 5 listing 7
Relative and absolute comparison routine. Classes for conveniently disabling and re-enabling (or vice- versa) floating-point
exceptions.
bool AlmostEqualRelativeAndAbs(float A, float B,
float maxDiff, float maxRelDiff = FLT_EPSILON) // Declare an object of this type in a scope in order to enable a
{ // specified set of floating-point exceptions temporarily. The old
// Check if the numbers are really close -- needed // exception state will be reset at the end.
// when comparing numbers near zero. // This class can be nested.
float diff = fabs(A - B); class FPExceptionEnabler
if (diff <= maxDiff) {
return true; public:
// Overflow, divide-by-zero, and invalid-operation are the FP
A = fabs(A); // exceptions most frequently associated with bugs.
B = fabs(B); FPExceptionEnabler(unsigned int enableBits = _EM_OVERFLOW |
float largest = (B > A) ? B : A; _EM_ZERODIVIDE | _EM_INVALID)
{
if (diff <= largest * maxRelDiff) // Retrieve the current state of the exception flags. This
return true; // must be done before changing them. _MCW_EM is a bit
return false; // mask representing all available exception masks.
} _controlfp_s(&mOldValues, _MCW_EM, _MCW_EM);
listing 6 // Make sure no non-exception flags have been specified,
// to avoid accidental changing of rounding modes, etc.
Relative and absolute comparison routine using ULPs. enableBits &= _MCW_EM;
bool AlmostEqualUlpsAndAbs(float A, float B, // Clear any pending FP exceptions. This must be done
float maxDiff, int maxUlpsDiff = 1) // prior to enabling FP exceptions since otherwise there
{ // may be a “deferred crash” as soon the exceptions are
// Check if the numbers are really close -- needed // enabled.
// when comparing numbers near zero. _clearfp();
float absDiff = fabs(A - B);
if (absDiff <= maxDiff) // Zero out the specified bits, leaving other bits alone.
return true; _controlfp_s(0, ~enableBits, enableBits);
}
Float_t uA(A); ~FPExceptionEnabler()
Float_t uB(B); {
// Reset the exception state.
// Different signs means they do not match. _controlfp_s(0, mOldValues, _MCW_EM);
if (uA.Negative() != uB.Negative()) }
return false;
private:
// Find the difference in ULPs. unsigned int mOldValues;
int ulpsDiff = abs(uA.i - uB.i);
if (ulpsDiff <= maxUlpsDiff) // Make the copy constructor and assignment operator private
return true; // and unimplemented to prohibit copying.
FPExceptionEnabler(const FPExceptionEnabler&);
return false; FPExceptionEnabler& operator=(const FPExceptionEnabler&);
} };

calculation should be zero, but Most people are inclined for values of x near pi, sin(x) error in float(pi) (and hence
float math stubbornly refuses to blame the implementation gives us the distance between x in sin(float(pi)) is roughly pi *
to give this answer. Typical of the sin function, but that is and pi. In other words, when we FLT_EPSILON / 2, assuming that
results for sin(float(pi)) are inappropriate—the problem isn’t calculate sin(float(pi)), we are float(pi) has an error of no more
-0.000000087422776, which with the sin function (which is actually calculating pi – float(pi). than half of one ULP. If we use
is 867,941,678 ULPs away from extremely accurate) but with The result is an extremely this as the absolute epsilon, then
zero. That’s as many ULPs as the value that is passed to it. It accurate measure of the error in the result suddenly seems totally
1.5e31 is from 1.0. Which is is important to understand that float(pi). A bit of calculus shows reasonable, and is “equal” to
not close. float(pi) is not equal to pi, and us that the expected maximum zero. The table below shows how

38 gamE DEvElOPER | OCTObER 2012